CAMT--------------------------------------------------------------------
CAMT  DFT_D Calculates the Grimme Dispersion Correction for the 
CAMT  energy and molecular gradient. 
CAMT  
CAMT  The code is tested against DFT-D3 V3.0 Rev 0 by S. Grimme
CAMT
CAMT  DFT-D2
CAMT  ------
CAMT  S. Grimme, J. Comput. Chem., 27 (2006), 1787-1799
CAMT
CAMT  DFT-D3
CAMT  ------
CAMT  S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys. 132 (2010), 154104
CAMT
CAMT  DFT-D3 + BJ Damping
CAMT  -------------------
CAMT  S. Grimme, S. Ehrlich and L. Goerigk, J. Comput. Chem. 32 (2011), 1456-1465
CAMT 
CAMT  The DFT-D3 empirical parameters are take from the DFT-D3 code v3.0 rev 0.
CAMT
CAMT----------------------------------------------------------------------

      SUBROUTINE DFT_D_DAL_IFC(EDISP,NDERIV,WORK,LWORK)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
#include "pgroup.h"
#include "chrxyz.h"
#include "chrsgn.h"
#include "chrnos.h"

      DIMENSION WORK(LWORK)
CAMT----------------------------------------------------------------------
CAMT  Variable info:
CAMT  
CAMT  EDISP     - The dispersion energy correction (output)
CAMT  NDERIV    - The number of derivatives to calculate (presently 
CAMT              max 1 -- add hessian later?) (input)
CAMT  N_ELEM    - the element Z value upto which the D3 parameters are
CAMT              defined
CAMT  N_ELEM_D2 - the element Z value upto which the D2 parameters are
CAMT              defined 
CAMT  MX_CORD   - Maximum number of coordination number references per
CAMT              element 
CAMT
CAMT  NOTE:       The gradient if requested is placed in a common block
CAMT              with the other gradient contributions -- see energy.h
CAMT              GRADFTD(). (optional output via common)
CAMT----------------------------------------------------------------------
      N_ELEM=94
      N_ELEM_D2=86
      MX_CORD=5

CAMT Perform the memory allocations... (all in one go!)
      KC6AB = 1
      KR0AB = KC6AB + N_ELEM*N_ELEM*MX_CORD*MX_CORD*3
      KFCOORD = KR0AB + N_ELEM*N_ELEM
      KTEMPG = KFCOORD + NATOMS*3
      KFCOORD2 = KTEMPG + NATOMS*3
      KCHARGE2 = KFCOORD2 + NATOMS*3
      KMXC = KCHARGE2 + NATOMS
      KNBASE = KMXC + N_ELEM
      KNSIGN = KNBASE + 8
      KR0 = KNSIGN + 8
      KC6 = KR0 + N_ELEM_D2
      KRCOV = KC6 + N_ELEM_D2
      KR2R4 = KRCOV + N_ELEM 
      KICOMP = KR2R4 + N_ELEM
      KR2AB = KICOMP + NATOMS*NATOMS
      KCC6AB = KR2AB + NATOMS*NATOMS
      KDMP = KCC6AB + NATOMS*NATOMS
      KSKIPDFD = KDMP + NATOMS*NATOMS
      KDRIJ = KSKIPDFD + (NATOMS*(NATOMS+1)/2) 
      KDC6_SAV = KDRIJ + (NATOMS*(NATOMS+1)/2)
      KDC6IJ = KDC6_SAV + (NATOMS*(NATOMS+1)/2)
      KDCN = KDC6IJ + NATOMS*NATOMS
      KD2 = KDCN + (NATOMS*(NATOMS+1)/2)
      KLAST = KD2 +3 

      IF (KLAST.GT.LWORK) CALL QUIT('ERROR: Insufficient Mem. in DFT_D')
      LWRK = LWORK-KLAST

      CALL DFT_D(EDISP,NDERIV,
     &           WORK(KC6AB),WORK(KR0AB),WORK(KFCOORD),
     &           WORK(KTEMPG),WORK(KFCOORD2),WORK(KCHARGE2),
     &           WORK(KMXC),WORK(KNBASE),WORK(KNSIGN),
     &           WORK(KR0),WORK(KC6),WORK(KRCOV),WORK(KR2R4),
     &           WORK(KICOMP),WORK(KR2AB),
     &           WORK(KCC6AB),WORK(KDMP),WORK(KSKIPDFD),
     &           WORK(KDRIJ),WORK(KDC6_SAV),WORK(KDC6IJ),
     &           WORK(KDCN), WORK(KD2),
     &           N_ELEM,N_ELEM_D2,MX_CORD,
     &           WORK,LWRK) 

      RETURN
      END

      SUBROUTINE DFT_D(EDISP,NDERIV,C6AB, R0AB, FCOORD,
     &                 TEMPG, FCOORD2, CHARGE2, MXC, NBASE, NSIGN,
     &                 R0, C6, RCOV, R2R4, ICOMP, R2AB,
     &                 CC6AB, DMP, SKIPDFD, DRIJ, DC6_SAV, DC6IJ,
     &                 DCN,D2, N_ELEM, N_ELEM_D2, MX_CORD,
     &                 WORK,LWORK) 
#include "implicit.h" 
#include "priunit.h" 
#include "maxaqn.h" 
#include "mxcent.h" 
#include "maxorb.h" 
#include "cbinuc.h" 
#include "nuclei.h" 
#include "symmet.h" 
#include "energy.h" 
#include "dftcom.h" 
#include "dftd.h" 
#include "pgroup.h"
#include "chrxyz.h"
#include "chrsgn.h"
#include "chrnos.h"

      DIMENSION WORK(LWORK)
!C6AB has dimensions of No. elements x No. elements x Coord No. Refs x Coord No. Refs. x 3
      DIMENSION C6AB(N_ELEM,N_ELEM,MX_CORD,MX_CORD,3)      
      DIMENSION R0AB(N_ELEM,N_ELEM)            
      DIMENSION FCOORD(NATOMS*3)       
      DIMENSION TEMPG(NATOMS*3)        
      DIMENSION FCOORD2(NATOMS,3)      
      DIMENSION CHARGE2(NATOMS)        
      INTEGER MXC(N_ELEM)                  
      INTEGER NBASE(8),NSIGN(8)        
      DIMENSION R0(N_ELEM_D2),C6(N_ELEM_D2)         
      DIMENSION RCOV(N_ELEM)    
      DIMENSION R2R4(N_ELEM)     
      INTEGER ICOMP(NATOMS*NATOMS)  
      DOUBLE PRECISION R2AB(NATOMS*NATOMS), CC6AB(NATOMS*NATOMS),
     &                 DMP(NATOMS*NATOMS) 
      DOUBLE PRECISION D2(3) 
      LOGICAL SKIPDFD(NATOMS*(NATOMS+1)/2) 
      DOUBLE PRECISION DRIJ(NATOMS*(NATOMS+1)/2)       
      DOUBLE PRECISION DC6_SAV(NATOMS*(NATOMS+1)/2)    
      DOUBLE PRECISION DCN(NATOMS*(NATOMS+1)/2)        
      DOUBLE PRECISION DC6IJ(NATOMS,NATOMS)            

      logical L_GRAD
CAMT External c-code routines for functional dependent parameter setup
      !DFT-D2
      external DFT_D2_S6, DFT_D2_RS6, DFT_D2_ALP 
      !DFT-D3
      external DFT_D3_S6, DFT_D3_RS6, DFT_D3_ALP, DFT_D3_RS18, 
     &         DFT_D3_S18
      !DFT-D3-BJ
      external DFT_D3BJ_S6, DFT_D3BJ_RS6, DFT_D3BJ_ALP, DFT_D3BJ_RS18, 
     &         DFT_D3BJ_S18

CAMT Parameters k1, k2 and k3
      double precision k1, k2, k3
      parameter (k1=16.0d0)
      parameter (k2=4.0d0/3.0d0)
      parameter (k3=-4.0d0)

CAMT Get the coordinates and charges, in particular if using symmetry, unpack from
CAMT symmetry coordinates

      CALL GET_COORD_CHAR(FCOORD,FCOORD2,CHARGE2)

CAMT  Determine which flavour of DFT-D we are to run (2,3,3+BJ)
      IF (DO_DFTD2) THEN
        NVER=2
      ELSE IF (DO_DFTD3) THEN
        NVER=3
      ENDIF

      CALL HEADER('DFT-D Empirical Dispersion Correction',-1)
      WRITE(LUPRI,'(A31,I5)')'Running DFT-D Version : ',NVER
      IF (DO_DFTD2) THEN
       WRITE(LUPRI,'(A31,A49)')'Reference : ',
     & 'S. Grimme, J. Comput. Chem., 27 (2006), 1787-1799'
      ELSE IF (DO_DFTD3) THEN
       WRITE(LUPRI,'(A31,A47)')'Reference : ',
     & 'S. Grimme, J. Antony, S. Ehrlich and H. Krieg, '
       WRITE(LUPRI,'(A31,A33)')' ',
     & 'J. Chem. Phys. 132 (2010), 154104'
        IF (DO_BJDAMP) THEN
        WRITE(LUPRI,'(A31,A38)')'Reference : ',
     & 'S. Grimme, S. Ehrlich and L. Goerigk, '
        WRITE(LUPRI,'(A31,A37)')' ',
     & 'J. Comput. Chem. 32 (2011), 1456-1465'
        ENDIF
      ENDIF

      IF (DO_DFTD3) THEN
        WRITE(LUPRI,'(A31,3X,L1)')'Becke-Johnson Damping : ',
     &                                                         DO_BJDAMP
        WRITE(LUPRI,'(A31,3X,L1)')'3-body terms : ',
     &                                                         DO_3BODY
      ENDIF

CAMT  Check not Becke-Johnson damping / 3-body terms when running with DFT-D2
      IF ((DO_DFTD2 .AND. DO_BJDAMP).OR.(DO_DFTD2 .AND. DO_3BODY)) THEN
       CALL QUIT('ERROR: INPUT PROBLEM IN DFT_D')
      ENDIF

CAMT  Determine what functional parameters are relevant and their values
      IF (DO_DFTD2) THEN
         IF (L_INP_D2PAR) THEN
           WRITE(LUPRI,*)' '
           WRITE(LUPRI,'(A55)')'** User Input DFT-D2 Parameters **'
           WRITE(LUPRI,*)' '
           S6  = D2_s6_inp
           RS6 = D2_rs6_inp
           ALP = D2_alp_inp
           S18 = 0.0d0
         ELSE
           S6  = DFT_D2_S6()
           RS6 = DFT_D2_RS6()
           ALP = DFT_D2_ALP()
           S18 = 0.0d0
         ENDIF
         CALL HEADER('DFT-D2 Functional Dependent Parameters',-1)
         WRITE(LUPRI,'(2(3X,A8,F20.12))')'s_6 = ',S6,'sr_6 = ',RS6
         WRITE(LUPRI,'(2(3X,A8,F20.12))')'s_8 = ',S18,'alp = ',alp
      ENDIF

      IF (DO_DFTD3) THEN
        IF (L_INP_D3PAR) THEN
         WRITE(LUPRI,*)' '
         WRITE(LUPRI,'(A55)')'** User Input DFT-D3 Parameters **'
         WRITE(LUPRI,*)' '
         S6   = D3_s6_inp
         ALP  = D3_alp_inp
         RS6  = D3_rs6_inp
         S18  = D3_s18_inp
         RS18 = D3_rs18_inp
        ELSE
         IF (DO_BJDAMP) THEN
           S6   = DFT_D3BJ_S6() 
           ALP  = DFT_D3BJ_ALP()
           RS6  = DFT_D3BJ_RS6()
           S18  = DFT_D3BJ_S18()
           RS18 = DFT_D3BJ_RS18()
         ELSE
           S6   = DFT_D3_S6()
           ALP  = DFT_D3_ALP()
           RS6  = DFT_D3_RS6()
           S18  = DFT_D3_S18()
           RS18 = DFT_D3_RS18()
         ENDIF
        ENDIF
        CALL HEADER('DFT-D3 Functional Dependent Parameters',-1)
        WRITE(LUPRI,'(2(3X,A9,F20.12))')'s_6 = ',S6,'sr_6 = ',RS6
        WRITE(LUPRI,'(2(3X,A9,F20.12))')'s_8 = ',S18,'sr_8 = ',RS18
        WRITE(LUPRI,'(3X,A9,F20.12)')'alp = ',alp
      ENDIF

CAMT  Set the run type (currently only energy / gradients)
      IF (NDERIV.EQ.0) THEN
        L_GRAD=.FALSE.
      ELSE IF (NDERIV.EQ.1) THEN
        L_GRAD=.TRUE.
        IF (L_GRAD.AND.DO_3BODY) THEN
           CALL QUIT('ERROR: Grad. for 3-body Corrections NYI, E Only')
        ENDIF
      ELSE
        CALL QUIT('ERROR: DFT_D - E and G only, use 1st Ord. Opts.')
      ENDIF

CAMT  Calculate the Energy Correction

      IF (DO_DFTD2) THEN
        CALL DFTD2_ENERGY(FCOORD2,CHARGE2,EDISP,S6,ALP,
     &                    RS6,
     &                 N_ELEM,N_ELEM_D2,MX_CORD,R0,C6)
      ENDIF

      IF (DO_DFTD3) THEN
        ALP6 = ALP
        ALP8 = ALP+2.0d0
        CALL DFTD3_ENERGY(FCOORD2,CHARGE2,S6,S18,EDISP,RS6,
     &                    RS18,ALP6,ALP8,E6,E8,C6AB,
     &                    R0AB,MXC,
     &                    ICOMP,R2AB,CC6AB,DMP,D2,
     &                    N_ELEM,N_ELEM_D2,MX_CORD,
     &                    RCOV,R2R4)
      ENDIF
 
      CALL AROUND('Dispersion Energy Correction')
      WRITE(LUPRI,'(A31,3X,F20.12)')'E_disp : ',EDISP
      
CAMT  Calculate the Gradient Contributions
      IF (L_GRAD) THEN
        IF (DO_DFTD2) THEN
          CALL DFTD2_GRAD(FCOORD2,CHARGE2,S6,ALP,RS6,
     &                 N_ELEM,N_ELEM_D2,MX_CORD,R0,C6)
        ENDIF

        IF (DO_DFTD3) THEN
          IF (DO_BJDAMP) THEN
            CALL DFTD3_GRAD_BJ(FCOORD2,CHARGE2,S6,S18,RS6,RS18,
     &                         ALP6,ALP8,E6,E8,C6AB,
     &                         R0AB,MXC,
     &                         SKIPDFD,
     &                         DRIJ,DC6_SAV,DCN,DC6IJ,
     &                         N_ELEM,N_ELEM_D2,MX_CORD,
     &                         RCOV,R2R4)   
          ELSE
            CALL DFTD3_GRAD(FCOORD2,CHARGE2,S6,S18,RS6,RS18,
     &                      ALP6,ALP8,E6,E8,C6AB,
     &                      R0AB,MXC,
     &                      SKIPDFD,
     &                      DRIJ,DC6_SAV,DCN,DC6IJ,
     &                      N_ELEM,N_ELEM_D2,MX_CORD,
     &                      RCOV,R2R4)
          ENDIF
        ENDIF

        IPRLOC=0
        IF (IPRLOC.GT.1 .OR. DFTD_TEST) THEN
          CALL HEADER('Cart. DFT-D Gradient Contribution',-1)
          IOFF = 0
          DO IATOM = 1, NUCDEP
            WRITE (LUPRI, '(1X,A6,F17.10,2F24.10)')
     &             NAMDEP(IATOM), (GRADFTD(IOFF+J),J=1,3)
                   IOFF = IOFF + 3
          END DO
          IF (DFTD_TEST) THEN
              WRITE (LUPRI,*) 'DFT-D TEST: Completed'
              CALL QUIT_NO_ERROR('DFT-D TEST completed')
          ENDIF
        ENDIF

CAMT If we are running with symmetry re-pack the cartesian gradient
CAMT to symmetry coordinates
        IF (MAXREP.GT.0) THEN
          CALL PACK_SYM_GRAD(TEMPG,NSIGN,NBASE,WORK,LWORK)
        ENDIF
      ENDIF

      RETURN
      END


CAMT ------------------------------ DFT-D2 ENERGY -------------------------------------
      SUBROUTINE DFTD2_ENERGY(FCOORD2,CHARGE2,EDISP,S6,ALP,
     &                        RS6,N_ELEM,N_ELEM_D2,MX_CORD,
     &                        R0,C6)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
      DIMENSION R0(N_ELEM_D2),C6(N_ELEM_D2)         
      DIMENSION FCOORD2(NATOMS,3)     
      DIMENSION CHARGE2(NATOMS)       
C     convert angstrom to bohr 
      PARAMETER (ang_to_au=1.88972612d0)  
C     convert Joule*nm^6/mol to Bohr^6*hartree 
C        1 Bohr = 52.9177 * 10^-3 nm 
C        1 Hartree = 2.6255*10^6 Joule/mol 
      PARAMETER (c6_conv=17.3452771d0) 
C     Initialize
      E = 0.0D0 
C
C     Setup parameters -- for now we assume runs are for E or (E+G) so these pass
C                         on to the gradient run. i.e. never just G
      CALL SET_R0(N_ELEM_D2,R0)
      CALL SET_C6(N_ELEM_D2,C6)
C
C     Loop over pairs of atoms A-B
C 
C     Run over nuclei A 
C 
      DO NCENTA = 1, NATOMS-1 
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM_D2 .AND. ICHARGA .GT. 0) THEN
            CORDAX = FCOORD2(NCENTA,1) 
            CORDAY = FCOORD2(NCENTA,2) 
            CORDAZ = FCOORD2(NCENTA,3) 
            C6A    = C6(ICHARGA)*c6_conv
            RvdWA  = R0(ICHARGA)*ang_to_au

C 
C           Run over nuclei B 
C 
            DO NCENTB =  NCENTA+1, NATOMS 
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM_D2 .AND. ICHARGB .GT. 0) THEN
                  CORDBX = FCOORD2(NCENTB,1) 
                  CORDBY = FCOORD2(NCENTB,2) 
                  CORDBZ = FCOORD2(NCENTB,3) 
                  C6B    = C6(ICHARGB)*c6_conv
                  RvdWB  = R0(ICHARGB)*ang_to_au
 
C                 C6 factor for atom pair A-B 
                  C6FAC = SQRT(C6A*C6B) 

C                 distance R between atoms A and B 
C                 and R^6  
                  RX    = CORDAX-CORDBX 
                  RY    = CORDAY-CORDBY 
                  RZ    = CORDAZ-CORDBZ 
                  R2    = RX*RX + RY*RY + RZ*RZ 
                  R6FAC = R2*R2*R2 
                  R     = SQRT(R2) 
 
C                 sum of the van der Waals radii RR   
                  RR   = RS6*(RvdWA + RvdWB) 
 
C                 damping function FDMP 
                  ALPHA = -ALP*R/RR+ALP 
                  EXPOA = EXP(ALPHA) 
                  FDMP  = 1.0D0/(1.0D0 + EXPOA) 
 
C                 dispersion correction contribution 
                  EADD  = C6FAC/R6FAC * FDMP 
                  E = E + EADD 
               END IF 
            END DO   
         END IF 
      END DO 
C     final dispersion correction 
C     S6 factor is functional dependent and passed in 
      EDISP = -S6 * E 
      RETURN 
      END                                                              

CAMT ------------------------------ DFT-D2 GRAD -------------------------------------
      SUBROUTINE DFTD2_GRAD(FCOORD2,CHARGE2,S6,ALP,RS6,
     &                      N_ELEM,N_ELEM_D2,MX_CORD,
     &                      R0,C6)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
      DIMENSION R0(N_ELEM_D2),C6(N_ELEM_D2)       
      DIMENSION FCOORD2(NATOMS,3)   
      DIMENSION CHARGE2(NATOMS)     
C     convert angstrom to bohr 
      PARAMETER (ang_to_au=1.88972612d0)
C     convert Joule*nm^6/mol to Bohr^6*hartree 
C        1 Bohr = 52.9177 * 10^-3 nm 
C        1 Hartree = 2.6255*10^6 Joule/mol 
      PARAMETER (c6_conv=17.3452771d0)

      R2_THR=9000.0d0   
      R2_THR2=1600.0d0

C     Initialize
      E = 0.0D0
      CALL DZERO(GRADFTD,MXCOOR)
C
C     Loop over pairs of atoms A-B
C 
C     Run over nuclei A 
C 
      DO NCENTA = 1, NATOMS-1
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM_D2 .AND. ICHARGA .GT. 0) THEN
            ISCOOA = (NCENTA-1)*3 
            CORDAX = FCOORD2(NCENTA,1)
            CORDAY = FCOORD2(NCENTA,2)
            CORDAZ = FCOORD2(NCENTA,3)
            C6A    = C6(ICHARGA)*c6_conv
            RvdWA  = R0(ICHARGA)*ang_to_au

C 
C           Run over nuclei B 
C 
            DO NCENTB =  NCENTA+1, NATOMS 
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM_D2 .AND. ICHARGB .GT. 0) THEN
                  ISCOOB = (NCENTB-1)*3 
                  CORDBX = FCOORD2(NCENTB,1)
                  CORDBY = FCOORD2(NCENTB,2)
                  CORDBZ = FCOORD2(NCENTB,3)
                  C6B    = C6(ICHARGB)*c6_conv
                  RvdWB  = R0(ICHARGB)*ang_to_au

C                 C6 factor for atom pair A-B 
                  C6FAC = SQRT(C6A*C6B)*S6

C                 distance R between atoms A and B 
C                 and R^6  
                  RX    = CORDAX-CORDBX
                  RY    = CORDAY-CORDBY
                  RZ    = CORDAZ-CORDBZ
                  R2    = RX*RX + RY*RY + RZ*RZ
                  R6FAC = R2*R2*R2
                  R     = SQRT(R2)
                  R7FAC = R6FAC*R

C                 sum of the van der Waals radii RR   
                  RR   = RS6*(RvdWA + RvdWB)

C                 damping function FDMP 
                  ALPHA = -ALP*R/RR+ALP
                  EXPOA = EXP(ALPHA) !EXPOA=DAMP6

                  FDMP  = (1.0D0 + EXPOA) !FDMP=DAMP1

                  T1 = EXPOA/(FDMP*FDMP*R7FAC*RR)
                  T2 = 6.0d0/(FDMP*R*R7FAC)
C  Gradient Contribs from first atom in pair
                  GRDX1 = (ALP*RX*T1-T2*RX)
                  GRDY1 = (ALP*RY*T1-T2*RY)
                  GRDZ1 = (ALP*RZ*T1-T2*RZ)
C  Gradient Contribs from second atom in pair
                  GRDX2 = (ALP*(-RX)*T1+T2*RX)
                  GRDY2 = (ALP*(-RY)*T1+T2*RY)
                  GRDZ2 = (ALP*(-RZ)*T1+T2*RZ)
C  GRADFTD is in energy.h (dimension MXCOOR)
                  GRADFTD(ISCOOA+1) = GRADFTD(ISCOOA+1) - GRDX1*C6FAC 
                  GRADFTD(ISCOOA+2) = GRADFTD(ISCOOA+2) - GRDY1*C6FAC 
                  GRADFTD(ISCOOA+3) = GRADFTD(ISCOOA+3) - GRDZ1*C6FAC 
                  GRADFTD(ISCOOB+1) = GRADFTD(ISCOOB+1) - GRDX2*C6FAC 
                  GRADFTD(ISCOOB+2) = GRADFTD(ISCOOB+2) - GRDY2*C6FAC 
                  GRADFTD(ISCOOB+3) = GRADFTD(ISCOOB+3) - GRDZ2*C6FAC 

               END IF
            END DO  
         END IF
      END DO
      RETURN
      END

CAMT ------------------------------ DFT-D3 ENERGY -------------------------------------
      SUBROUTINE DFTD3_ENERGY(FCOORD2,CHARGE2,S6,S8,EDISP,RS6,RS8,ALP6,
     &                        ALP8,E6,E8,C6AB,R0AB,MXC,
     &                        ICOMP,R2AB,CC6AB,DMP,D2,
     &                        N_ELEM,N_ELEM_D2,MX_CORD,
     &                        RCOV,R2R4)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
      DIMENSION RCOV(N_ELEM)     
      DIMENSION R2R4(N_ELEM)     
      DIMENSION C6AB(N_ELEM,N_ELEM,MX_CORD,MX_CORD,3) 
      DIMENSION R0AB(N_ELEM,N_ELEM)       
      DIMENSION FCOORD2(NATOMS,3) 
      DIMENSION CHARGE2(NATOMS)   

      INTEGER ICOMP(NATOMS*NATOMS)  
      DOUBLE PRECISION R2AB(NATOMS*NATOMS), CC6AB(NATOMS*NATOMS),  
     &                 DMP(NATOMS*NATOMS) 
      DOUBLE PRECISION D2(3) 
      INTEGER MXC(N_ELEM) 
      INTEGER GET_INDX
C     convert angstrom to bohr 
      PARAMETER (ang_to_au=1.88972612d0)  
C     convert Joule*nm^6/mol to Bohr^6*hartree 
C        1 Bohr = 52.9177 * 10^-3 nm 
C        1 Hartree = 2.6255*10^6 Joule/mol 
      PARAMETER (c6_conv=17.3452771d0) 

      R2_THR=9000.0d0   
      R2_THR2=1600.0d0
      CN_THRESH = 1600.0d0

      IPRLOC=0
      IF (IPRLOC.GT.1) THEN
        CALL HEADER('Cutoff Parameters',-1)
        WRITE(LUPRI,'(A30,F20.12)')'Cut-off:',DSQRT(R2_THR)
        WRITE(LUPRI,'(A30,F20.12)')'Coord No. Cutoff:',DSQRT(CN_THRESH)
      ENDIF

C     Initialize
      E6 = 0.0D0 
      E8 = 0.0D0 
      E63 = 0.0D0
C
C Setup the parameters  -- Here we assume runs are for E or E+G, 
C                          so paras are set here and passed onto
C                          the gradient run 
      CALL SET_RCOV(N_ELEM,RCOV)
      CALL SET_R2R4(N_ELEM,R2R4)
      CALL SET_R0AB(N_ELEM,ang_to_au,R0AB)
      CALL LOAD_C6AB(MX_CORD,N_ELEM,c6ab,mxc)

C     Set parameters for Becke-Johnson Damping
      bj_par1 = RS6
      bj_par2 = RS8

C     Calculate the coordinaton numbers (set in CN in CBLOCK)
      CALL SET_COORD_NOS(FCOORD2,CHARGE2,RCOV,N_ELEM)

C     Loop over pairs of atoms A-B
C 
C     Run over nuclei A 
C 
      DO NCENTA = 1, NATOMS-1 
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM .AND. ICHARGA .GT. 0) THEN
            CORDAX = FCOORD2(NCENTA,1) 
            CORDAY = FCOORD2(NCENTA,2) 
            CORDAZ = FCOORD2(NCENTA,3) 
C 
C           Run over nuclei B 
C 
            DO NCENTB =  NCENTA+1, NATOMS 
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM .AND. ICHARGB .GT. 0) THEN
                  CORDBX = FCOORD2(NCENTB,1) 
                  CORDBY = FCOORD2(NCENTB,2) 
                  CORDBZ = FCOORD2(NCENTB,3) 

C                 distance R between atoms A and B 
C                 and R^6  
                  RX    = CORDAX-CORDBX 
                  RY    = CORDAY-CORDBY 
                  RZ    = CORDAZ-CORDBZ 
                  R2    = RX*RX + RY*RY + RZ*RZ 

                  IF (R2.LT.R2_THR) THEN
                    R      = SQRT(R2) 
                    R6     = R2**3 
                    R8     = R6*R2 
 
C                   sum of the van der Waals radii RR   
 
                    RR = R0AB(ICHARGB,ICHARGA) / R
 
C                   damping function FDMP 
C                   c6 part
                    FDMP6  = 1.0D0/(1.0D0 + 6.0d0*(RS6*RR)**ALP6) 
C                   c8 part
                    FDMP8  = 1.0D0/(1.0D0 + 6.0d0*(RS8*RR)**ALP8) 

C                   get c6 coeffient
                    CALL SET_C6_COEFF(CHARGE2,MXC,C6AB,NCENTA,NCENTB,C6,
     &                                N_ELEM,MX_CORD)
C                   calculate c8 coefficient
                    C8 = 3.0d0*C6*R2R4(ICHARGA)*R2R4(ICHARGB)

C                   if requested to add 3-body terms store the required parts
                    IF (DO_3BODY) THEN
                      IJ = GET_INDX(NCENTB,NCENTA)
                      ICOMP(IJ) = 1
                      CC6AB(IJ) = SQRT(C6)
                      R2AB(IJ)  = R2
                      DMP(IJ) = (1.0D0/RR)**(1.0d0/3.0d0)
                    END IF                   

C                   accumulate the E6 and E8 dispersion contributions 
                    IF (DO_BJDAMP) THEN
                      SRC8C6 = SQRT(C8/C6)
                      E6  = E6 + C6 / (R6+(bj_par1*SRC8C6+bj_par2)**6) 
                      E8  = E8 + C8 / (R8+(bj_par1*SRC8C6+bj_par2)**8) 
                    ELSE
                      E6  = E6 + C6*FDMP6/R6
                      E8  = E8 + C8*FDMP8/R8 
                    ENDIF
                  END IF 
               END IF 
            END DO   
         END IF 
      END DO

C     Construct 3-body corrections and add
      IF (DO_3BODY) THEN
C 
C     Run over nuclei A 
C 
      DO NCENTA = 1, NATOMS-1
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM .AND. ICHARGA .GT. 0) THEN
            CORDAX = FCOORD2(NCENTA,1)
            CORDAY = FCOORD2(NCENTA,2)
            CORDAZ = FCOORD2(NCENTA,3)
C 
C           Run over nuclei B 
C 
            DO NCENTB =  NCENTA+1, NATOMS 
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM .AND. ICHARGB .GT. 0) THEN
                  CORDBX = FCOORD2(NCENTB,1)
                  CORDBY = FCOORD2(NCENTB,2)
                  CORDBZ = FCOORD2(NCENTB,3)

                  IJ = GET_INDX(NCENTB,NCENTA)

                  IF (ICOMP(IJ).EQ.1) THEN
                    DO NCENTC = NCENTB+1, NATOMS 
                      IK=GET_INDX(NCENTC,NCENTA)
                      JK=GET_INDX(NCENTC,NCENTB)
                      IF ((ICOMP(IK).EQ.1).AND.(ICOMP(JK).EQ.1)) THEN
                        !damping function product
                        RAV = (4.0D0/3.0D0)/(DMP(IK)*DMP(JK)*DMP(IJ))
                        FDMP = 1.D0/( 1.D0+6.D0*RAV**ALP8 )
                        !triple C6 coefficient
                        C9 = CC6AB(IJ)*CC6AB(IK)*CC6AB(JK)
                        !angular terms
                        d2(1) = r2ab(ij)
                        d2(2) = r2ab(jk)
                        d2(3) = r2ab(ik)
                        t1 = (d2(1)+d2(2)-d2(3))/sqrt(d2(1)*d2(2))
                        t2 = (d2(1)+d2(3)-d2(2))/sqrt(d2(1)*d2(3))
                        t3 = (d2(3)+d2(2)-d2(1))/sqrt(d2(2)*d2(3))
                        ang=0.375d0*t1*t2*t3+1.0d0
                        !compute E63
                      E63 = E63 - FDMP*C9*ANG/(d2(1)*d2(2)*d2(3))**1.5D0
                      ENDIF
                    ENDDO
                  END IF
               END IF
            END DO
         END IF
      END DO

      ENDIF

      E6 = -E6*S6
      E8 = -E8*S8

      E63 = E63*S6

      EDISP = E6+E8+E63
 
      RETURN
      END                                                               

      INTEGER FUNCTION GET_INDX(I,J)
      INTEGER I,J,IMAX,IMIN
      IMAX=max(I,J)
      IMIN=min(I,J)
      GET_INDX=IMIN+IMAX*(IMAX-1)/2
      RETURN
      END


CAMT ------------------------------ DFT-D3 GRAD -------------------------------------
      SUBROUTINE DFTD3_GRAD(FCOORD2,CHARGE2,S6,S18,RS6,RS8,ALP6,
     &                      ALP8,E6,E8,C6AB,R0AB,MXC,
     &                      SKIPDFD,
     &                      DRIJ,DC6_SAV,DCN,DC6IJ,
     &                      N_ELEM,N_ELEM_D2,MX_CORD,
     &                      RCOV,R2R4)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
      LOGICAL SKIPDFD(NATOMS*(NATOMS+1)/2) !To allocate
      DIMENSION C6AB(N_ELEM,N_ELEM,MX_CORD,MX_CORD,3)   
      DIMENSION R0AB(N_ELEM,N_ELEM)         
      DIMENSION RCOV(N_ELEM),R2R4(N_ELEM)    
      INTEGER MXC(N_ELEM)               

      DIMENSION FCOORD2(NATOMS,3)    
      DIMENSION CHARGE2(NATOMS)     

      DOUBLE PRECISION DRIJ(NATOMS*(NATOMS+1)/2)       
      DOUBLE PRECISION DC6_SAV(NATOMS*(NATOMS+1)/2)    
      DOUBLE PRECISION DCN(NATOMS*(NATOMS+1)/2)        
      DOUBLE PRECISION DC6IJ(NATOMS,NATOMS)            
      INTEGER GET_INDX
      PARAMETER (ang_to_au=1.88972612d0)  

C  Parameters k1, k2 and k3
      double precision k1, k2, k3
      parameter (k1=16.0d0)
      parameter (k2=4.0d0/3.0d0)
      parameter (k3=-4.0d0)

      R2_THR=9000.0d0   
      R2_THR2=1600.0d0
      CN_THRESH = 1600.0d0

C  Initialize
      CALL DZERO(GRADFTD,MXCOOR)
      S8  = S18
      S10 = S18

      DISP=0
      DO I = 1, NATOMS*(NATOMS+1)/2
        DRIJ(I)=0.0d0
        DC6_SAV(I)=0.0d0
        DCN(I)=0.0d0
      ENDDO

      DO I = 1, NATOMS*(NATOMS+1)/2
        SKIPDFD(I) = .TRUE.
      ENDDO
C 
C     Run over nuclei A 
C 
      DO NCENTA = 1, NATOMS 
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM .AND. ICHARGA .GT. 0) THEN
            CORDAX = FCOORD2(NCENTA,1)
            CORDAY = FCOORD2(NCENTA,2)
            CORDAZ = FCOORD2(NCENTA,3)
C 
C           Run over nuclei B 
C 
            DO NCENTB =  1, NCENTA-1
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM .AND. ICHARGB .GT. 0) THEN
                  CORDBX = FCOORD2(NCENTB,1)
                  CORDBY = FCOORD2(NCENTB,2)
                  CORDBZ = FCOORD2(NCENTB,3)
C
C                 distance R between atoms A and B 
C                 and R^6  
                  RX    = CORDAX-CORDBX
                  RY    = CORDAY-CORDBY
                  RZ    = CORDAZ-CORDBZ
                  R2    = RX*RX + RY*RY + RZ*RZ
                  R     = DSQRT(R2)
                  R6    = R2*R2*R2
                  R7    = R6*R
                  R8    = R6*R2
                  R9    = R8*R

C                 sum of the van der Waals radii RR   

                  RR     = R0AB(ICHARGB,ICHARGA) 
                  R42    = R2R4(ICHARGA)*R2R4(ICHARGB)
                  RCOVIJ = RCOV(ICHARGA) + RCOV(ICHARGB) 

                  IF (R2.LT.R2_THR) THEN
                  IJ = GET_INDX(NCENTB,NCENTA)
                  SKIPDFD(IJ) = .FALSE.

C Calculate dC6(i,j)/dCN(i) and dC6(i,j)/dCN(j)
C results in dC6ij
                  CALL GET_dC6_dCNij(MX_CORD,N_ELEM,C6AB,
     &                               MXC(ICHARGA),
     &                               MXC(ICHARGB), CN(NCENTA),
     &                               CN(NCENTB), ICHARGA, 
     &                               ICHARGB, NCENTA, NCENTB,
     &                               C6, DC6IJ(NCENTA,NCENTB),
     &                               DC6IJ(NCENTB,NCENTA) )
              
C Calculate the damping function derivatives
                  T6 = (R/(RS6*RR))**(-ALP6)
                  FDMP6 = 1.0D0 / (1.0D0+6.0D0*T6)
                  T8 = (R/(RS8*RR))**(-ALP8)
                  FDMP8 = 1.0D0 / (1.0D0+6.0D0*T8)

                  ! d(r^(-6)) / d(r_ij) ...
                  DRIJ(IJ) = DRIJ(IJ) - S6*(6.0D0/R7*C6*FDMP6)
     &                                - S8*(24.0D0/R9*C6*R42*FDMP8)
                  ! d(f_dmp)/d(r_ij) ...
                  DRIJ(IJ) = DRIJ(IJ)
     &                     + S6*C6/R7*6.0D0*ALP6*T6*FDMP6*FDMP6
     &                     + S8*C6*R42/R9*18.0D0*ALP8*T8*FDMP8*FDMP8

C Store part of the dC6 term for later
                  DC6_SAV(IJ) = S6/R6*FDMP6+3.0D0*S8*R42/R8*FDMP8
                  
C Calculate dCN(i)/dr_ij = dCN(j)/dr_ij
                  IF (R2.LT.CN_THRESH) THEN
                    EXPT = exp(-k1*(RCOVIJ/R-1.0d0))
                    DCN(IJ) = -k1*RCOVIJ*EXPT/
     &                         (R*R*(EXPT+1.0D0)*(EXPT+1.0D0)) 

C Calculate dC6/dCN * dCN/dr_ij = dC6/dr_ij
                    DC6 = (DC6IJ(NCENTA,NCENTB)+DC6IJ(NCENTB,NCENTA))
     &                  * DCN(IJ)
                    ! d(C6(IJ))/d(r_ij) 
                    DRIJ(IJ) = DRIJ(IJ) + DC6_SAV(IJ)*DC6
                  ELSE
                    DC6 = 0.0D0
                    SKIPDFD(IJ) = .TRUE.
                  ENDIF


C  Loop to calculate dC6(i,k)/dr_ij, dC6(j,k)/dr_ij, dC6(ij)/dr_ik, dC6(ij)/dr_jk
                  DO NCENTC = 1, NCENTB-1

                    IK = GET_INDX(NCENTA,NCENTC)
                    JK = GET_INDX(NCENTB,NCENTC)

                    IF (.NOT. SKIPDFD(IJ)) THEN
                      DRIJ(IJ) =DRIJ(IJ) +
     &                          DC6_SAV(IK)*DC6IJ(NCENTA,NCENTC)*DCN(IJ)
                      DRIJ(IJ) =DRIJ(IJ) +
     &                          DC6_SAV(JK)*DC6IJ(NCENTB,NCENTC)*DCN(IJ)
                    ENDIF

                    IF (.NOT. SKIPDFD(JK)) THEN
                      DRIJ(JK) =DRIJ(JK) +
     &                          DC6_SAV(IK)*DC6IJ(NCENTC,NCENTA)*DCN(JK)
                      DRIJ(JK) =DRIJ(JK) +
     &                          DC6_SAV(IJ)*DC6IJ(NCENTB,NCENTA)*DCN(JK)
                    ENDIF

                    IF (.NOT. SKIPDFD(IK)) THEN
                      DRIJ(IK) =DRIJ(IK) +
     &                          DC6_SAV(JK)*DC6IJ(NCENTC,NCENTB)*DCN(IK)
                      DRIJ(IK) =DRIJ(IK) +
     &                          DC6_SAV(IJ)*DC6IJ(NCENTA,NCENTB)*DCN(IK)
                    ENDIF


                  ENDDO  !NCENTC

                  ENDIF

               ENDIF
            ENDDO !NCENTB
         ENDIF
      ENDDO !NCENTA


C All dE/dr_ij calculated. Now calculate dE/dr_ij * dr_ij/dxyz_i

      DO NCENTA = 2, NATOMS 
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM .AND. ICHARGA .GT. 0) THEN
            ISCOOA = (NCENTA-1)*3
            CORDAX = FCOORD2(NCENTA,1)
            CORDAY = FCOORD2(NCENTA,2)
            CORDAZ = FCOORD2(NCENTA,3)

            DO NCENTB = 1, NCENTA-1
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM .AND. ICHARGB .GT. 0) THEN
                  ISCOOB = (NCENTB-1)*3 
                  CORDBX = FCOORD2(NCENTB,1)
                  CORDBY = FCOORD2(NCENTB,2)
                  CORDBZ = FCOORD2(NCENTB,3)

!                 distance R between atoms A and B 
!                 and R^6  
                  RX    = CORDAX-CORDBX
                  RY    = CORDAY-CORDBY
                  RZ    = CORDAZ-CORDBZ
                  R2    = RX*RX + RY*RY + RZ*RZ
                  R     = DSQRT(R2)

                  IJ = GET_INDX(NCENTA,NCENTB)
                  GRDX1 = -DRIJ(IJ)*RX/R
                  GRDY1 = -DRIJ(IJ)*RY/R
                  GRDZ1 = -DRIJ(IJ)*RZ/R
                  GRDX2 =  DRIJ(IJ)*RX/R
                  GRDY2 =  DRIJ(IJ)*RY/R
                  GRDZ2 =  DRIJ(IJ)*RZ/R

C GRADFTD is in energy.h (dimension MXCOOR)
                  GRADFTD(ISCOOA+1) = GRADFTD(ISCOOA+1) + GRDX1
                  GRADFTD(ISCOOA+2) = GRADFTD(ISCOOA+2) + GRDY1
                  GRADFTD(ISCOOA+3) = GRADFTD(ISCOOA+3) + GRDZ1
                  GRADFTD(ISCOOB+1) = GRADFTD(ISCOOB+1) + GRDX2
                  GRADFTD(ISCOOB+2) = GRADFTD(ISCOOB+2) + GRDY2
                  GRADFTD(ISCOOB+3) = GRADFTD(ISCOOB+3) + GRDZ2
               ENDIF          
            ENDDO
         ENDIF
      ENDDO

      RETURN
      END


CAMT ------------------------------ DFT-D3 GRAD with Becke-Johnson damping-------------------------------------
      SUBROUTINE DFTD3_GRAD_BJ(FCOORD2,CHARGE2,S6,S18,RS6,RS8,ALP6,
     &                         ALP8,E6,E8,C6AB,R0AB,MXC,
     &                         SKIPDFD,
     &                         DRIJ,DC6_SAV,DCN,DC6IJ,
     &                         N_ELEM,N_ELEM_D2,MX_CORD,
     &                         RCOV,R2R4)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
      LOGICAL SKIPDFD(NATOMS*(NATOMS+1)/2)    
      DIMENSION C6AB(N_ELEM,N_ELEM,MX_CORD,MX_CORD,3)             
      DIMENSION R0AB(N_ELEM,N_ELEM)                   
      DIMENSION RCOV(N_ELEM),R2R4(N_ELEM)             
      INTEGER MXC(N_ELEM)                         

      DIMENSION FCOORD2(NATOMS,3)             
      DIMENSION CHARGE2(NATOMS)               

      DOUBLE PRECISION DRIJ(NATOMS*(NATOMS+1)/2)      
      DOUBLE PRECISION DC6_SAV(NATOMS*(NATOMS+1)/2)   
      DOUBLE PRECISION DCN(NATOMS*(NATOMS+1)/2)       
      DOUBLE PRECISION DC6IJ(NATOMS,NATOMS)           
      INTEGER GET_INDX
      PARAMETER (ang_to_au=1.88972612d0)

C Parameters k1, k2 and k3
      double precision k1, k2, k3
      parameter (k1=16.0d0)
      parameter (k2=4.0d0/3.0d0)
      parameter (k3=-4.0d0)

      R2_THR=9000.0d0   
      R2_THR2=1600.0d0
      CN_THRESH = 1600.0d0

C Initialize
      CALL DZERO(GRADFTD,MXCOOR)
      A1 = RS6
      A2 = RS8
      S8  = S18

      DISP=0
      DO I = 1, NATOMS*(NATOMS+1)/2
        DRIJ(I)=0.0d0
        DC6_SAV(I)=0.0d0
        DCN(I)=0.0d0
      ENDDO

      DO I = 1, NATOMS*(NATOMS+1)/2
        SKIPDFD(I) = .TRUE.
      ENDDO
C 
C     Run over nuclei A 
C 
      DO NCENTA = 1, NATOMS 
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM .AND. ICHARGA .GT. 0) THEN
            CORDAX = FCOORD2(NCENTA,1)
            CORDAY = FCOORD2(NCENTA,2)
            CORDAZ = FCOORD2(NCENTA,3)
C 
C           Run over nuclei B 
C 
            DO NCENTB =  1, NCENTA-1
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM .AND. ICHARGB .GT. 0) THEN
                  CORDBX = FCOORD2(NCENTB,1)
                  CORDBY = FCOORD2(NCENTB,2)
                  CORDBZ = FCOORD2(NCENTB,3)

C                 distance R between atoms A and B 
C                 and R^6  
                  RX    = CORDAX-CORDBX
                  RY    = CORDAY-CORDBY
                  RZ    = CORDAZ-CORDBZ
                  R2    = RX*RX + RY*RY + RZ*RZ
                  R     = DSQRT(R2)
                  R4    = R2*R2
                  R6    = R2*R2*R2
                  R7    = R6*R
                  R8    = R6*R2
                  R9    = R8*R

C                 sum of the van der Waals radii RR   

                  RR     = R0AB(ICHARGB,ICHARGA)
                  R42    = R2R4(ICHARGA)*R2R4(ICHARGB)
                  RCOVIJ = RCOV(ICHARGA) + RCOV(ICHARGB)

                  IF (R2.LT.R2_THR) THEN
                  IJ = GET_INDX(NCENTB,NCENTA)
                  SKIPDFD(IJ) = .FALSE.

C Calculate dC6(i,j)/dCN(i) and dC6(i,j)/dCN(j)
C results in dC6ij
                  CALL GET_dC6_dCNij(MX_CORD,N_ELEM,C6AB,
     &                               MXC(ICHARGA),
     &                               MXC(ICHARGB), CN(NCENTA),
     &                               CN(NCENTB), ICHARGA,
     &                               ICHARGB, NCENTA, NCENTB,
     &                               C6, DC6IJ(NCENTA,NCENTB),
     &                               DC6IJ(NCENTB,NCENTA) )
C Calculate the damping functions and their derivatives
C Use the Becke--Johnson radius
                  R0BJ = A1*SQRT(3.0D0*R42)+A2   

                  T6 = (R6+R0BJ**6)
                  T8 = (R8+R0BJ**8)

                  DRIJ(IJ) = DRIJ(IJ)
     &                     - S6*C6*6.0D0*R4*R/(T6*T6)
     &                     - S8*C6*24.0D0*R42*R6*R/(T8*T8)

C Store part of the dC6 term for later
                  DC6_SAV(IJ) = S6/T6+3.0D0*S8*R42/T8

C Calculate dCN(i)/dr_ij = dCN(j)/dr_ij
                  IF (R2.LT.CN_THRESH) THEN
                    EXPT = exp(-k1*(RCOVIJ/R-1.0d0))
                    DCN(IJ) = -k1*RCOVIJ*EXPT/
     &                         (R*R*(EXPT+1.0D0)*(EXPT+1.0D0))

C Calculate dC6/dCN * dCN/dr_ij = dC6/dr_ij
                    DC6 = (DC6IJ(NCENTA,NCENTB)+DC6IJ(NCENTB,NCENTA))
     &                  * DCN(IJ)
                    ! d(C6(IJ))/d(r_ij) 
                    DRIJ(IJ) = DRIJ(IJ) + DC6_SAV(IJ)*DC6
                  ELSE
                    DC6 = 0.0D0
                    DCN(IJ) = 0.0d0
                    SKIPDFD(IJ) = .TRUE.
                  ENDIF

C Loop to calculate dC6(i,k)/dr_ij, dC6(j,k)/dr_ij, dC6(ij)/dr_ik, dC6(ij)/dr_jk
                  DO NCENTC = 1, NCENTB-1

                    IK = GET_INDX(NCENTA,NCENTC)
                    JK = GET_INDX(NCENTB,NCENTC)


                    DC6=DC6IJ(NCENTA,NCENTC)*DCN(IJ)
                    DRIJ(IJ) = DRIJ(IJ) + DC6_SAV(IK)*DC6

                    DC6 = DC6IJ(NCENTB,NCENTC)*DCN(IJ)
                    DRIJ(IJ) = DRIJ(IJ) + DC6_SAV(JK)*DC6

                    IF (.NOT. SKIPDFD(JK)) THEN
                      DRIJ(JK) =DRIJ(JK) +
     &                          DC6_SAV(IK)*DC6IJ(NCENTC,NCENTA)*DCN(JK)
                      DRIJ(JK) =DRIJ(JK) +
     &                          DC6_SAV(IJ)*DC6IJ(NCENTB,NCENTA)*DCN(JK)
                    ENDIF

                    IF (.NOT. SKIPDFD(IK)) THEN
                      DRIJ(IK) =DRIJ(IK) +
     &                          DC6_SAV(JK)*DC6IJ(NCENTC,NCENTB)*DCN(IK)
                      DRIJ(IK) =DRIJ(IK) +
     &                          DC6_SAV(IJ)*DC6IJ(NCENTA,NCENTB)*DCN(IK)
                    ENDIF


                  ENDDO  !NCENTC

                  ENDIF

               ENDIF
            ENDDO !NCENTB
         ENDIF
      ENDDO !NCENTA

C All dE/dr_ij calculated. Now calculate dE/dr_ij * dr_ij/dxyz_i

      DO NCENTA = 2, NATOMS 
         ICHARGA = NINT( CHARGE2(NCENTA) )
         IF (ICHARGA .LE. N_ELEM .AND. ICHARGA .GT. 0) THEN
            ISCOOA = (NCENTA-1)*3
            CORDAX = FCOORD2(NCENTA,1)
            CORDAY = FCOORD2(NCENTA,2)
            CORDAZ = FCOORD2(NCENTA,3)

            DO NCENTB = 1, NCENTA-1
               ICHARGB = NINT( CHARGE2(NCENTB) )
               IF (ICHARGB .LE. N_ELEM .AND. ICHARGB .GT. 0) THEN
                  ISCOOB = (NCENTB-1)*3
                  CORDBX = FCOORD2(NCENTB,1)
                  CORDBY = FCOORD2(NCENTB,2)
                  CORDBZ = FCOORD2(NCENTB,3)

C                 distance R between atoms A and B 
C                 and R^6  
                  RX    = CORDAX-CORDBX
                  RY    = CORDAY-CORDBY
                  RZ    = CORDAZ-CORDBZ
                  R2    = RX*RX + RY*RY + RZ*RZ
                  R     = DSQRT(R2)

                  IJ = GET_INDX(NCENTA,NCENTB)
                  GRDX1 = -DRIJ(IJ)*RX/R
                  GRDY1 = -DRIJ(IJ)*RY/R
                  GRDZ1 = -DRIJ(IJ)*RZ/R
                  GRDX2 =  DRIJ(IJ)*RX/R
                  GRDY2 =  DRIJ(IJ)*RY/R
                  GRDZ2 =  DRIJ(IJ)*RZ/R

C GRADFTD is in energy.h (dimension MXCOOR)
                  GRADFTD(ISCOOA+1) = GRADFTD(ISCOOA+1) + GRDX1
                  GRADFTD(ISCOOA+2) = GRADFTD(ISCOOA+2) + GRDY1
                  GRADFTD(ISCOOA+3) = GRADFTD(ISCOOA+3) + GRDZ1
                  GRADFTD(ISCOOB+1) = GRADFTD(ISCOOB+1) + GRDX2
                  GRADFTD(ISCOOB+2) = GRADFTD(ISCOOB+2) + GRDY2
                  GRADFTD(ISCOOB+3) = GRADFTD(ISCOOB+3) + GRDZ2
               ENDIF
            ENDDO
         ENDIF
      ENDDO

      RETURN
      END

!---------------------- DFT-D3 SET COORDINATION NUMBERS --------------------------------

      SUBROUTINE SET_COORD_NOS(FCOORD2,CHARGE2,RCOV,N_ELEM)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
      dimension RCOV(N_ELEM)
      DIMENSION FCOORD2(NATOMS,3)
      DIMENSION CHARGE2(NATOMS)
C Parameters k1, k2 and k3
      double precision k1, k2, k3
      parameter (k1=16.0d0)
      parameter (k2=4.0d0/3.0d0)
      parameter (k3=-4.0d0)
      CN_THRESH = 1600.0d0

      DO I = 1,NATOMS
        CN_I=0.0d0
        DO J = 1, NATOMS
          IF (I.NE.J) THEN
            RX = FCOORD2(J,1)-FCOORD2(I,1)
            RY = FCOORD2(J,2)-FCOORD2(I,2)
            RZ = FCOORD2(J,3)-FCOORD2(I,3)
            R2    = RX*RX + RY*RY + RZ*RZ
            IF (R2.LT.CN_THRESH) THEN
C Covalent separation in au
              R  = SQRT(R2)
              IZ=INT(CHARGE2(I))
              JZ=INT(CHARGE2(J))
              RRCOV = RCOV(IZ) + RCOV(JZ)
              RR = RRCOV/R
C Coord-No Damping Function
              CNDAMP = 1.0d0 / (1.0d0+exp(-k1*(RR-1.0d0)))
              CN_I = CN_I + CNDAMP
            ENDIF
          ENDIF
        ENDDO
        CN(I)=CN_I
      ENDDO
      RETURN
      END


!-------------------- DFT-D3 Set C6 parameters -----------------------------------------
      SUBROUTINE SET_C6_COEFF(CHARGE2,MXC,C6AB,NCENTA,NCENTB,C6,
     &                        N_ELEM,MX_CORD)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
      DOUBLE PRECISION C6AB(N_ELEM,N_ELEM,MX_CORD,MX_CORD,3)
      DIMENSION CHARGE2(NATOMS)
      INTEGER MXC(N_ELEM)
C Parameters k1, k2 and k3
      double precision k1, k2, k3
      parameter (k1=16.0d0)
      parameter (k2=4.0d0/3.0d0)
      parameter (k3=-4.0d0)
      parameter (BIGN = 1.0D99)
      integer   CHRGA, CHRGB

      C6MEM  = -BIGN
      RSUM   = 0.0D0
      CSUM   = 0.0D0
      C6     = 0.0D0
      R_SAVE = BIGN

      CHRGA=INT(CHARGE2(NCENTA))
      CHRGB=INT(CHARGE2(NCENTB))

C Coord No.s in common block 
      CN_I=CN(NCENTA)
      CN_J=CN(NCENTB)

      DO I = 1, MXC(CHRGA)            !Loop over coord no. refs of atom 1
        DO J = 1, MXC(CHRGB)          !Loop over coord no. refs of atom 2
          C6 = C6AB(CHRGA,CHRGB,I,J,1)
          IF (C6.GT.0) THEN
            CN1 = C6AB(CHRGA,CHRGB,I,J,2) 
            CN2 = C6AB(CHRGA,CHRGB,I,J,3)
            !distance
            R = (CN1-CN_I)**2 + (CN2-CN_J)**2
            IF (R.LT.R_SAVE) THEN
               R_SAVE = R
               C6MEM  = C6
            ENDIF
            TMP1 = exp(k3*R)
            RSUM = RSUM + TMP1
            CSUM = CSUM + TMP1*C6
          ENDIF
        ENDDO
      ENDDO

      IF (RSUM.GT.1.0D-99) THEN
         C6 = CSUM/RSUM
      ELSE
         C6 = C6MEM
      END IF

      RETURN
      END

      SUBROUTINE LOAD_C6AB(MX_CORD,N_ELEM,C6AB,MXCi)
      IMPLICIT NONE
      INTEGER MX_CORD,N_ELEM,MXCi(N_ELEM),NLINES
      REAL*8  C6AB(N_ELEM,N_ELEM,MX_CORD,MX_CORD,3)
      INTEGER IAT,JAT,L,J,K,IL,IADR,JADR,NN,KK
#include "dft_dftd_pars.h"
      !nlines set in the above file
      C6AB=-1
      MXCi=0
      KK=1
      DO NN=1,NLINES
       IAT=INT(PARS(KK+1))
       JAT=INT(PARS(KK+2))
       CALL LIMIT(IAT,JAT,IADR,JADR)
       MXCi(IAT)=MAX(MXCi(IAT),IADR)
       MXCi(JAT)=MAX(MXCi(JAT),JADR)

       C6AB(IAT,JAT,IADR,JADR,1)=PARS(KK)
       C6AB(IAT,JAT,IADR,JADR,2)=PARS(KK+3)
       C6AB(IAT,JAT,IADR,JADR,3)=PARS(KK+4)

       C6AB(JAT,IAT,JADR,IADR,1)=PARS(KK)
       C6AB(JAT,IAT,JADR,IADR,2)=PARS(KK+4)
       C6AB(JAT,IAT,JADR,IADR,3)=PARS(KK+3)
       KK=(NN*5)+1
      ENDDO
      RETURN
      END  


      SUBROUTINE LIMIT(IAT,JAT,IADR,JADR)
      IMPLICIT NONE
      INTEGER IAT,JAT,IADR,JADR,I
      IADR=1
      JADR=1
      I=100
 10   IF(IAT.GT.100) THEN
         IAT=IAT-100
         IADR=IADR+1
         GOTO 10
      ENDIf

      I=100
 20   IF(JAT.GT.100) THEN
         JAT=JAT-100
         JADR=JADR+1
         GOTO 20
      ENDIf
      RETURN
      END 


      SUBROUTINE GET_dC6_dCNij(MX_CORD,N_ELEM,C6AB,MXCI,MXCJ,CNI,CNJ,
     &           IZI,IZJ,IAT,JAT,C6CHECK,DC6I,DC6J)

      IMPLICIT NONE
! Parameters k1, k2 and k3
      DOUBLE PRECISION K1, K2, K3
      PARAMETER (K1=16.0d0)
      PARAMETER (K2=4.0d0/3.0d0)
      PARAMETER (K3=-4.0d0)
      INTEGER MX_CORD,N_ELEM
      REAL*8 C6AB(N_ELEM,N_ELEM,MX_CORD,MX_CORD,3)
      INTEGER MXCI,MXCJ   !mxc(iz(iat))
      REAL*8 CNI,CNJ
      INTEGER IAT,JAT,IZI,IZJ
      REAL*8  DC6I,DC6J,C6CHECK

      INTEGER I,J,A,B
      REAL*8 NUMER,DENOM,DNUMER_I,DDENOM_I,DNUMER_J,DDENOM_J
      REAL*8 EXPTERM,CN_REFI,CN_REFJ,C6REF,R
      REAL*8 C6MEM,R_SAVE

      C6MEM=-1.d99
      R_SAVE=9999.0
      NUMER=0.0d0
      DENOM=0.0d0

      DNUMER_I=0.d0
      DDENOM_I=0.d0
      DNUMER_J=0.d0
      DDENOM_J=0.d0


      DO A=1,MXCI
        DO B=1,MXCJ
          C6REF=C6AB(IZI,IZJ,A,B,1)
          IF (C6REF.GT.0) THEN
            CN_REFI=C6AB(IZI,IZJ,A,B,2)
            CN_REFJ=C6AB(IZI,IZJ,A,B,3)
            R=(CN_REFI-CNI)*(CN_REFI-CNI)+(CN_REFJ-CNJ)*(CN_REFJ-CNJ)
            IF (R.LT.R_SAVE) THEN
               R_SAVE=R
               C6MEM=C6REF
            ENDIF
            EXPTERM=EXP(K3*R)
            NUMER=NUMER+C6REF*EXPTERM
            DENOM=DENOM+EXPTERM
            DNUMER_I=DNUMER_I+C6REF*EXPTERM*
     &             2.D0*K3*(CNI-CN_REFI)
            DDENOM_I=DDENOM_I+EXPTERM*
     &             2.D0*K3*(CNI-CN_REFI)

            DNUMER_J=DNUMER_J+C6REF*EXPTERM*
     &             2.D0*K3*(CNJ-CN_REFJ)
            DDENOM_J=DDENOM_J+EXPTERM*
     &             2.D0*K3*(CNJ-CN_REFJ)
          ENDIF
        ENDDO !b
      ENDDO !a

      IF (DENOM.GT.1.0D-99) THEN
        C6CHECK=NUMER/DENOM
        DC6I=((DNUMER_I*DENOM)-(DDENOM_I*NUMER))
     &    /(DENOM*DENOM)
        DC6J=((DNUMER_J*DENOM)-(DDENOM_J*NUMER))
     &    /(DENOM*DENOM)
      ELSE
        C6CHECK=C6MEM
        DC6I=0.0D0
        DC6J=0.0D0
      ENDIF
      END


      SUBROUTINE SET_R0AB(N_ELEM,ANG_TO_AU,R)
      IMPLICIT NONE
      INTEGER N_ELEM,I,J,K
      REAL*8 R(N_ELEM,N_ELEM),ANG_TO_AU
      REAL*8 R0AB(4465)
      R0AB(   1:  70)=(/
     &   2.1823,  1.8547,  1.7347,  2.9086,  2.5732,  3.4956,  2.3550
     &,  2.5095,  2.9802,  3.0982,  2.5141,  2.3917,  2.9977,  2.9484
     &,  3.2160,  2.4492,  2.2527,  3.1933,  3.0214,  2.9531,  2.9103
     &,  2.3667,  2.1328,  2.8784,  2.7660,  2.7776,  2.7063,  2.6225
     &,  2.1768,  2.0625,  2.6395,  2.6648,  2.6482,  2.5697,  2.4846
     &,  2.4817,  2.0646,  1.9891,  2.5086,  2.6908,  2.6233,  2.4770
     &,  2.3885,  2.3511,  2.2996,  1.9892,  1.9251,  2.4190,  2.5473
     &,  2.4994,  2.4091,  2.3176,  2.2571,  2.1946,  2.1374,  2.9898
     &,  2.6397,  3.6031,  3.1219,  3.7620,  3.2485,  2.9357,  2.7093
     &,  2.5781,  2.4839,  3.7082,  2.5129,  2.7321,  3.1052,  3.2962
     &/)
      R0AB(  71: 140)=(/
     &   3.1331,  3.2000,  2.9586,  3.0822,  2.8582,  2.7120,  3.2570
     &,  3.4839,  2.8766,  2.7427,  3.2776,  3.2363,  3.5929,  3.2826
     &,  3.0911,  2.9369,  2.9030,  2.7789,  3.3921,  3.3970,  4.0106
     &,  2.8884,  2.6605,  3.7513,  3.1613,  3.3605,  3.3325,  3.0991
     &,  2.9297,  2.8674,  2.7571,  3.8129,  3.3266,  3.7105,  3.7917
     &,  2.8304,  2.5538,  3.3932,  3.1193,  3.1866,  3.1245,  3.0465
     &,  2.8727,  2.7664,  2.6926,  3.4608,  3.2984,  3.5142,  3.5418
     &,  3.5017,  2.6190,  2.4797,  3.1331,  3.0540,  3.0651,  2.9879
     &,  2.9054,  2.8805,  2.7330,  2.6331,  3.2096,  3.5668,  3.3684
     &,  3.3686,  3.3180,  3.3107,  2.4757,  2.4019,  2.9789,  3.1468
     &/)
      R0AB( 141: 210)=(/
     &   2.9768,  2.8848,  2.7952,  2.7457,  2.6881,  2.5728,  3.0574
     &,  3.3264,  3.3562,  3.2529,  3.1916,  3.1523,  3.1046,  2.3725
     &,  2.3289,  2.8760,  2.9804,  2.9093,  2.8040,  2.7071,  2.6386
     &,  2.5720,  2.5139,  2.9517,  3.1606,  3.2085,  3.1692,  3.0982
     &,  3.0352,  2.9730,  2.9148,  3.2147,  2.8315,  3.8724,  3.4621
     &,  3.8823,  3.3760,  3.0746,  2.8817,  2.7552,  2.6605,  3.9740
     &,  3.6192,  3.6569,  3.9586,  3.6188,  3.3917,  3.2479,  3.1434
     &,  4.2411,  2.7597,  3.0588,  3.3474,  3.6214,  3.4353,  3.4729
     &,  3.2487,  3.3200,  3.0914,  2.9403,  3.4972,  3.7993,  3.6773
     &,  3.8678,  3.5808,  3.8243,  3.5826,  3.4156,  3.8765,  4.1035
     &/)
      R0AB( 211: 280)=(/
     &   2.7361,  2.9765,  3.2475,  3.5004,  3.4185,  3.4378,  3.2084
     &,  3.2787,  3.0604,  2.9187,  3.4037,  3.6759,  3.6586,  3.8327
     &,  3.5372,  3.7665,  3.5310,  3.3700,  3.7788,  3.9804,  3.8903
     &,  2.6832,  2.9060,  3.2613,  3.4359,  3.3538,  3.3860,  3.1550
     &,  3.2300,  3.0133,  2.8736,  3.4024,  3.6142,  3.5979,  3.5295
     &,  3.4834,  3.7140,  3.4782,  3.3170,  3.7434,  3.9623,  3.8181
     &,  3.7642,  2.6379,  2.8494,  3.1840,  3.4225,  3.2771,  3.3401
     &,  3.1072,  3.1885,  2.9714,  2.8319,  3.3315,  3.5979,  3.5256
     &,  3.4980,  3.4376,  3.6714,  3.4346,  3.2723,  3.6859,  3.8985
     &,  3.7918,  3.7372,  3.7211,  2.9230,  2.6223,  3.4161,  2.8999
     &/)
      R0AB( 281: 350)=(/
     &   3.0557,  3.3308,  3.0555,  2.8508,  2.7385,  2.6640,  3.5263
     &,  3.0277,  3.2990,  3.7721,  3.5017,  3.2751,  3.1368,  3.0435
     &,  3.7873,  3.2858,  3.2140,  3.1727,  3.2178,  3.4414,  2.5490
     &,  2.7623,  3.0991,  3.3252,  3.1836,  3.2428,  3.0259,  3.1225
     &,  2.9032,  2.7621,  3.2490,  3.5110,  3.4429,  3.3845,  3.3574
     &,  3.6045,  3.3658,  3.2013,  3.6110,  3.8241,  3.7090,  3.6496
     &,  3.6333,  3.0896,  3.5462,  2.4926,  2.7136,  3.0693,  3.2699
     &,  3.1272,  3.1893,  2.9658,  3.0972,  2.8778,  2.7358,  3.2206
     &,  3.4566,  3.3896,  3.3257,  3.2946,  3.5693,  3.3312,  3.1670
     &,  3.5805,  3.7711,  3.6536,  3.5927,  3.5775,  3.0411,  3.4885
     &/)
      R0AB( 351: 420)=(/
     &   3.4421,  2.4667,  2.6709,  3.0575,  3.2357,  3.0908,  3.1537
     &,  2.9235,  3.0669,  2.8476,  2.7054,  3.2064,  3.4519,  3.3593
     &,  3.2921,  3.2577,  3.2161,  3.2982,  3.1339,  3.5606,  3.7582
     &,  3.6432,  3.5833,  3.5691,  3.0161,  3.4812,  3.4339,  3.4327
     &,  2.4515,  2.6338,  3.0511,  3.2229,  3.0630,  3.1265,  2.8909
     &,  3.0253,  2.8184,  2.6764,  3.1968,  3.4114,  3.3492,  3.2691
     &,  3.2320,  3.1786,  3.2680,  3.1036,  3.5453,  3.7259,  3.6090
     &,  3.5473,  3.5327,  3.0018,  3.4413,  3.3907,  3.3593,  3.3462
     &,  2.4413,  2.6006,  3.0540,  3.1987,  3.0490,  3.1058,  2.8643
     &,  2.9948,  2.7908,  2.6491,  3.1950,  3.3922,  3.3316,  3.2585
     &/)
      R0AB( 421: 490)=(/
     &   3.2136,  3.1516,  3.2364,  3.0752,  3.5368,  3.7117,  3.5941
     &,  3.5313,  3.5164,  2.9962,  3.4225,  3.3699,  3.3370,  3.3234
     &,  3.3008,  2.4318,  2.5729,  3.0416,  3.1639,  3.0196,  3.0843
     &,  2.8413,  2.7436,  2.7608,  2.6271,  3.1811,  3.3591,  3.3045
     &,  3.2349,  3.1942,  3.1291,  3.2111,  3.0534,  3.5189,  3.6809
     &,  3.5635,  3.5001,  3.4854,  2.9857,  3.3897,  3.3363,  3.3027
     &,  3.2890,  3.2655,  3.2309,  2.8502,  2.6934,  3.2467,  3.1921
     &,  3.5663,  3.2541,  3.0571,  2.9048,  2.8657,  2.7438,  3.3547
     &,  3.3510,  3.9837,  3.6871,  3.4862,  3.3389,  3.2413,  3.1708
     &,  3.6096,  3.6280,  3.6860,  3.5568,  3.4836,  3.2868,  3.3994
     &/)
      R0AB( 491: 560)=(/
     &   3.3476,  3.3170,  3.2950,  3.2874,  3.2606,  3.9579,  2.9226
     &,  2.6838,  3.7867,  3.1732,  3.3872,  3.3643,  3.1267,  2.9541
     &,  2.8505,  2.7781,  3.8475,  3.3336,  3.7359,  3.8266,  3.5733
     &,  3.3959,  3.2775,  3.1915,  3.9878,  3.8816,  3.5810,  3.5364
     &,  3.5060,  3.8097,  3.3925,  3.3348,  3.3019,  3.2796,  3.2662
     &,  3.2464,  3.7136,  3.8619,  2.9140,  2.6271,  3.4771,  3.1774
     &,  3.2560,  3.1970,  3.1207,  2.9406,  2.8322,  2.7571,  3.5455
     &,  3.3514,  3.5837,  3.6177,  3.5816,  3.3902,  3.2604,  3.1652
     &,  3.7037,  3.6283,  3.5858,  3.5330,  3.4884,  3.5789,  3.4094
     &,  3.3473,  3.3118,  3.2876,  3.2707,  3.2521,  3.5570,  3.6496
     &/)
      r0ab( 561: 630)=(/
     &   3.6625,  2.7300,  2.5870,  3.2471,  3.1487,  3.1667,  3.0914
     &,  3.0107,  2.9812,  2.8300,  2.7284,  3.3259,  3.3182,  3.4707
     &,  3.4748,  3.4279,  3.4182,  3.2547,  3.1353,  3.5116,  3.9432
     &,  3.8828,  3.8303,  3.7880,  3.3760,  3.7218,  3.3408,  3.3059
     &,  3.2698,  3.2446,  3.2229,  3.4422,  3.5023,  3.5009,  3.5268
     &,  2.6026,  2.5355,  3.1129,  3.2863,  3.1029,  3.0108,  2.9227
     &,  2.8694,  2.8109,  2.6929,  3.1958,  3.4670,  3.4018,  3.3805
     &,  3.3218,  3.2815,  3.2346,  3.0994,  3.3937,  3.7266,  3.6697
     &,  3.6164,  3.5730,  3.2522,  3.5051,  3.4686,  3.4355,  3.4084
     &,  3.3748,  3.3496,  3.3692,  3.4052,  3.3910,  3.3849,  3.3662
     &/)
      R0AB( 631: 700)=(/
     &   2.5087,  2.4814,  3.0239,  3.1312,  3.0535,  2.9457,  2.8496
     &,  2.7780,  2.7828,  2.6532,  3.1063,  3.3143,  3.3549,  3.3120
     &,  3.2421,  3.1787,  3.1176,  3.0613,  3.3082,  3.5755,  3.5222
     &,  3.4678,  3.4231,  3.1684,  3.3528,  3.3162,  3.2827,  3.2527
     &,  3.2308,  3.2029,  3.3173,  3.3343,  3.3092,  3.2795,  3.2452
     &,  3.2096,  3.2893,  2.8991,  4.0388,  3.6100,  3.9388,  3.4475
     &,  3.1590,  2.9812,  2.8586,  2.7683,  4.1428,  3.7911,  3.8225
     &,  4.0372,  3.7059,  3.4935,  3.3529,  3.2492,  4.4352,  4.0826
     &,  3.9733,  3.9254,  3.8646,  3.9315,  3.7837,  3.7465,  3.7211
     &,  3.7012,  3.6893,  3.6676,  3.7736,  4.0660,  3.7926,  3.6158
     &/)
      R0AB( 701: 770)=(/
     &   3.5017,  3.4166,  4.6176,  2.8786,  3.1658,  3.5823,  3.7689
     &,  3.5762,  3.5789,  3.3552,  3.4004,  3.1722,  3.0212,  3.7241
     &,  3.9604,  3.8500,  3.9844,  3.7035,  3.9161,  3.6751,  3.5075
     &,  4.1151,  4.2877,  4.1579,  4.1247,  4.0617,  3.4874,  3.9848
     &,  3.9280,  3.9079,  3.8751,  3.8604,  3.8277,  3.8002,  3.9981
     &,  3.7544,  4.0371,  3.8225,  3.6718,  4.3092,  4.4764,  2.8997
     &,  3.0953,  3.4524,  3.6107,  3.6062,  3.5783,  3.3463,  3.3855
     &,  3.1746,  3.0381,  3.6019,  3.7938,  3.8697,  3.9781,  3.6877
     &,  3.8736,  3.6451,  3.4890,  3.9858,  4.1179,  4.0430,  3.9563
     &,  3.9182,  3.4002,  3.8310,  3.7716,  3.7543,  3.7203,  3.7053
     &/)
      R0AB( 771: 840)=(/
     &   3.6742,  3.8318,  3.7631,  3.7392,  3.9892,  3.7832,  3.6406
     &,  4.1701,  4.3016,  4.2196,  2.8535,  3.0167,  3.3978,  3.5363
     &,  3.5393,  3.5301,  3.2960,  3.3352,  3.1287,  2.9967,  3.6659
     &,  3.7239,  3.8070,  3.7165,  3.6368,  3.8162,  3.5885,  3.4336
     &,  3.9829,  4.0529,  3.9584,  3.9025,  3.8607,  3.3673,  3.7658
     &,  3.7035,  3.6866,  3.6504,  3.6339,  3.6024,  3.7708,  3.7283
     &,  3.6896,  3.9315,  3.7250,  3.5819,  4.1457,  4.2280,  4.1130
     &,  4.0597,  3.0905,  2.7998,  3.6448,  3.0739,  3.2996,  3.5262
     &,  3.2559,  3.0518,  2.9394,  2.8658,  3.7514,  3.2295,  3.5643
     &,  3.7808,  3.6931,  3.4723,  3.3357,  3.2429,  4.0280,  3.5589
     &/)
      R0AB( 841: 910)=(/
     &   3.4636,  3.4994,  3.4309,  3.6177,  3.2946,  3.2376,  3.2050
     &,  3.1847,  3.1715,  3.1599,  3.5555,  3.8111,  3.7693,  3.5718
     &,  3.4498,  3.3662,  4.1608,  3.7417,  3.6536,  3.6154,  3.8596
     &,  3.0301,  2.7312,  3.5821,  3.0473,  3.2137,  3.4679,  3.1975
     &,  2.9969,  2.8847,  2.8110,  3.6931,  3.2076,  3.4943,  3.5956
     &,  3.6379,  3.4190,  3.2808,  3.1860,  3.9850,  3.5105,  3.4330
     &,  3.3797,  3.4155,  3.6033,  3.2737,  3.2145,  3.1807,  3.1596
     &,  3.1461,  3.1337,  3.4812,  3.6251,  3.7152,  3.5201,  3.3966
     &,  3.3107,  4.1128,  3.6899,  3.6082,  3.5604,  3.7834,  3.7543
     &,  2.9189,  2.6777,  3.4925,  2.9648,  3.1216,  3.2940,  3.0975
     &/)
      R0AB( 911: 980)=(/
     &   2.9757,  2.8493,  2.7638,  3.6085,  3.1214,  3.4006,  3.4793
     &,  3.5147,  3.3806,  3.2356,  3.1335,  3.9144,  3.4183,  3.3369
     &,  3.2803,  3.2679,  3.4871,  3.1714,  3.1521,  3.1101,  3.0843
     &,  3.0670,  3.0539,  3.3890,  3.5086,  3.5895,  3.4783,  3.3484
     &,  3.2559,  4.0422,  3.5967,  3.5113,  3.4576,  3.6594,  3.6313
     &,  3.5690,  2.8578,  2.6334,  3.4673,  2.9245,  3.0732,  3.2435
     &,  3.0338,  2.9462,  2.8143,  2.7240,  3.5832,  3.0789,  3.3617
     &,  3.4246,  3.4505,  3.3443,  3.1964,  3.0913,  3.8921,  3.3713
     &,  3.2873,  3.2281,  3.2165,  3.4386,  3.1164,  3.1220,  3.0761
     &,  3.0480,  3.0295,  3.0155,  3.3495,  3.4543,  3.5260,  3.4413
     &/)
      R0AB( 981:1050)=(/
     &   3.3085,  3.2134,  4.0170,  3.5464,  3.4587,  3.4006,  3.6027
     &,  3.5730,  3.4945,  3.4623,  2.8240,  2.5960,  3.4635,  2.9032
     &,  3.0431,  3.2115,  2.9892,  2.9148,  2.7801,  2.6873,  3.5776
     &,  3.0568,  3.3433,  3.3949,  3.4132,  3.3116,  3.1616,  3.0548
     &,  3.8859,  3.3719,  3.2917,  3.2345,  3.2274,  3.4171,  3.1293
     &,  3.0567,  3.0565,  3.0274,  3.0087,  2.9939,  3.3293,  3.4249
     &,  3.4902,  3.4091,  3.2744,  3.1776,  4.0078,  3.5374,  3.4537
     &,  3.3956,  3.5747,  3.5430,  3.4522,  3.4160,  3.3975,  2.8004
     &,  2.5621,  3.4617,  2.9154,  3.0203,  3.1875,  2.9548,  2.8038
     &,  2.7472,  2.6530,  3.5736,  3.0584,  3.3304,  3.3748,  3.3871
     &/)
      R0AB(1051:1120)=(/
     &   3.2028,  3.1296,  3.0214,  3.8796,  3.3337,  3.2492,  3.1883
     &,  3.1802,  3.4050,  3.0756,  3.0478,  3.0322,  3.0323,  3.0163
     &,  3.0019,  3.3145,  3.4050,  3.4656,  3.3021,  3.2433,  3.1453
     &,  3.9991,  3.5017,  3.4141,  3.3520,  3.5583,  3.5251,  3.4243
     &,  3.3851,  3.3662,  3.3525,  2.7846,  2.5324,  3.4652,  2.8759
     &,  3.0051,  3.1692,  2.9273,  2.7615,  2.7164,  2.6212,  3.5744
     &,  3.0275,  3.3249,  3.3627,  3.3686,  3.1669,  3.0584,  2.9915
     &,  3.8773,  3.3099,  3.2231,  3.1600,  3.1520,  3.4023,  3.0426
     &,  3.0099,  2.9920,  2.9809,  2.9800,  2.9646,  3.3068,  3.3930
     &,  3.4486,  3.2682,  3.1729,  3.1168,  3.9952,  3.4796,  3.3901
     &/)
      R0AB(1121:1190)=(/
     &   3.3255,  3.5530,  3.5183,  3.4097,  3.3683,  3.3492,  3.3360
     &,  3.3308,  2.5424,  2.6601,  3.2555,  3.2807,  3.1384,  3.1737
     &,  2.9397,  2.8429,  2.8492,  2.7225,  3.3875,  3.4910,  3.4520
     &,  3.3608,  3.3036,  3.2345,  3.2999,  3.1487,  3.7409,  3.8392
     &,  3.7148,  3.6439,  3.6182,  3.1753,  3.5210,  3.4639,  3.4265
     &,  3.4075,  3.3828,  3.3474,  3.4071,  3.3754,  3.3646,  3.3308
     &,  3.4393,  3.2993,  3.8768,  3.9891,  3.8310,  3.7483,  3.3417
     &,  3.3019,  3.2250,  3.1832,  3.1578,  3.1564,  3.1224,  3.4620
     &,  2.9743,  2.8058,  3.4830,  3.3474,  3.6863,  3.3617,  3.1608
     &,  3.0069,  2.9640,  2.8427,  3.5885,  3.5219,  4.1314,  3.8120
     &/)
      R0AB(1191:1260)=(/
     &   3.6015,  3.4502,  3.3498,  3.2777,  3.8635,  3.8232,  3.8486
     &,  3.7215,  3.6487,  3.4724,  3.5627,  3.5087,  3.4757,  3.4517
     &,  3.4423,  3.4139,  4.1028,  3.8388,  3.6745,  3.5562,  3.4806
     &,  3.4272,  4.0182,  3.9991,  4.0007,  3.9282,  3.7238,  3.6498
     &,  3.5605,  3.5211,  3.5009,  3.4859,  3.4785,  3.5621,  4.2623
     &,  3.0775,  2.8275,  4.0181,  3.3385,  3.5379,  3.5036,  3.2589
     &,  3.0804,  3.0094,  2.9003,  4.0869,  3.5088,  3.9105,  3.9833
     &,  3.7176,  3.5323,  3.4102,  3.3227,  4.2702,  4.0888,  3.7560
     &,  3.7687,  3.6681,  3.6405,  3.5569,  3.4990,  3.4659,  3.4433
     &,  3.4330,  3.4092,  3.8867,  4.0190,  3.7961,  3.6412,  3.5405
     &/)
      R0AB(1261:1330)=(/
     &   3.4681,  4.3538,  4.2136,  3.9381,  3.8912,  3.9681,  3.7909
     &,  3.6774,  3.6262,  3.5999,  3.5823,  3.5727,  3.5419,  4.0245
     &,  4.1874,  3.0893,  2.7917,  3.7262,  3.3518,  3.4241,  3.5433
     &,  3.2773,  3.0890,  2.9775,  2.9010,  3.8048,  3.5362,  3.7746
     &,  3.7911,  3.7511,  3.5495,  3.4149,  3.3177,  4.0129,  3.8370
     &,  3.7739,  3.7125,  3.7152,  3.7701,  3.5813,  3.5187,  3.4835
     &,  3.4595,  3.4439,  3.4242,  3.7476,  3.8239,  3.8346,  3.6627
     &,  3.5479,  3.4639,  4.1026,  3.9733,  3.9292,  3.8667,  3.9513
     &,  3.8959,  3.7698,  3.7089,  3.6765,  3.6548,  3.6409,  3.5398
     &,  3.8759,  3.9804,  4.0150,  2.9091,  2.7638,  3.5066,  3.3377
     &/)
      R0AB(1331:1400)=(/
     &   3.3481,  3.2633,  3.1810,  3.1428,  2.9872,  2.8837,  3.5929
     &,  3.5183,  3.6729,  3.6596,  3.6082,  3.5927,  3.4224,  3.2997
     &,  3.8190,  4.1865,  4.1114,  4.0540,  3.6325,  3.5697,  3.5561
     &,  3.5259,  3.4901,  3.4552,  3.4315,  3.4091,  3.6438,  3.6879
     &,  3.6832,  3.7043,  3.5557,  3.4466,  3.9203,  4.2919,  4.2196
     &,  4.1542,  3.7573,  3.7039,  3.6546,  3.6151,  3.5293,  3.4849
     &,  3.4552,  3.5192,  3.7673,  3.8359,  3.8525,  3.8901,  2.7806
     &,  2.7209,  3.3812,  3.4958,  3.2913,  3.1888,  3.0990,  3.0394
     &,  2.9789,  2.8582,  3.4716,  3.6883,  3.6105,  3.5704,  3.5059
     &,  3.4619,  3.4138,  3.2742,  3.7080,  3.9773,  3.9010,  3.8409
     &/)
      R0AB(1401:1470)=(/
     &   3.7944,  3.4465,  3.7235,  3.6808,  3.6453,  3.6168,  3.5844
     &,  3.5576,  3.5772,  3.5959,  3.5768,  3.5678,  3.5486,  3.4228
     &,  3.8107,  4.0866,  4.0169,  3.9476,  3.6358,  3.5800,  3.5260
     &,  3.4838,  3.4501,  3.4204,  3.3553,  3.6487,  3.6973,  3.7398
     &,  3.7405,  3.7459,  3.7380,  2.6848,  2.6740,  3.2925,  3.3386
     &,  3.2473,  3.1284,  3.0301,  2.9531,  2.9602,  2.8272,  3.3830
     &,  3.5358,  3.5672,  3.5049,  3.4284,  3.3621,  3.3001,  3.2451
     &,  3.6209,  3.8299,  3.7543,  3.6920,  3.6436,  3.3598,  3.5701
     &,  3.5266,  3.4904,  3.4590,  3.4364,  3.4077,  3.5287,  3.5280
     &,  3.4969,  3.4650,  3.4304,  3.3963,  3.7229,  3.9402,  3.8753
     &/)
      R0AB(1471:1540)=(/
     &   3.8035,  3.5499,  3.4913,  3.4319,  3.3873,  3.3520,  3.3209
     &,  3.2948,  3.5052,  3.6465,  3.6696,  3.6577,  3.6388,  3.6142
     &,  3.5889,  3.3968,  3.0122,  4.2241,  3.7887,  4.0049,  3.5384
     &,  3.2698,  3.1083,  2.9917,  2.9057,  4.3340,  3.9900,  4.6588
     &,  4.1278,  3.8125,  3.6189,  3.4851,  3.3859,  4.6531,  4.3134
     &,  4.2258,  4.1309,  4.0692,  4.0944,  3.9850,  3.9416,  3.9112
     &,  3.8873,  3.8736,  3.8473,  4.6027,  4.1538,  3.8994,  3.7419
     &,  3.6356,  3.5548,  4.8353,  4.5413,  4.3891,  4.3416,  4.3243
     &,  4.2753,  4.2053,  4.1790,  4.1685,  4.1585,  4.1536,  4.0579
     &,  4.1980,  4.4564,  4.2192,  4.0528,  3.9489,  3.8642,  5.0567
     &/)
      R0AB(1541:1610)=(/
     &   3.0630,  3.3271,  4.0432,  4.0046,  4.1555,  3.7426,  3.5130
     &,  3.5174,  3.2884,  3.1378,  4.1894,  4.2321,  4.1725,  4.1833
     &,  3.8929,  4.0544,  3.8118,  3.6414,  4.6373,  4.6268,  4.4750
     &,  4.4134,  4.3458,  3.8582,  4.2583,  4.1898,  4.1562,  4.1191
     &,  4.1069,  4.0639,  4.1257,  4.1974,  3.9532,  4.1794,  3.9660
     &,  3.8130,  4.8160,  4.8272,  4.6294,  4.5840,  4.0770,  4.0088
     &,  3.9103,  3.8536,  3.8324,  3.7995,  3.7826,  4.2294,  4.3380
     &,  4.4352,  4.1933,  4.4580,  4.2554,  4.1072,  5.0454,  5.1814
     &,  3.0632,  3.2662,  3.6432,  3.8088,  3.7910,  3.7381,  3.5093
     &,  3.5155,  3.3047,  3.1681,  3.7871,  3.9924,  4.0637,  4.1382
     &/)
      R0AB(1611:1680)=(/
     &   3.8591,  4.0164,  3.7878,  3.6316,  4.1741,  4.3166,  4.2395
     &,  4.1831,  4.1107,  3.5857,  4.0270,  3.9676,  3.9463,  3.9150
     &,  3.9021,  3.8708,  4.0240,  4.1551,  3.9108,  4.1337,  3.9289
     &,  3.7873,  4.3666,  4.5080,  4.4232,  4.3155,  3.8461,  3.8007
     &,  3.6991,  3.6447,  3.6308,  3.5959,  3.5749,  4.0359,  4.3124
     &,  4.3539,  4.1122,  4.3772,  4.1785,  4.0386,  4.7004,  4.8604
     &,  4.6261,  2.9455,  3.2470,  3.6108,  3.8522,  3.6625,  3.6598
     &,  3.4411,  3.4660,  3.2415,  3.0944,  3.7514,  4.0397,  3.9231
     &,  4.0561,  3.7860,  3.9845,  3.7454,  3.5802,  4.1366,  4.3581
     &,  4.2351,  4.2011,  4.1402,  3.5381,  4.0653,  4.0093,  3.9883
     &/)
      R0AB(1681:1750)=(/
     &   3.9570,  3.9429,  3.9112,  3.8728,  4.0682,  3.8351,  4.1054
     &,  3.8928,  3.7445,  4.3415,  4.5497,  4.3833,  4.3122,  3.8051
     &,  3.7583,  3.6622,  3.6108,  3.5971,  3.5628,  3.5408,  4.0780
     &,  4.0727,  4.2836,  4.0553,  4.3647,  4.1622,  4.0178,  4.5802
     &,  4.9125,  4.5861,  4.6201,  2.9244,  3.2241,  3.5848,  3.8293
     &,  3.6395,  3.6400,  3.4204,  3.4499,  3.2253,  3.0779,  3.7257
     &,  4.0170,  3.9003,  4.0372,  3.7653,  3.9672,  3.7283,  3.5630
     &,  4.1092,  4.3347,  4.2117,  4.1793,  4.1179,  3.5139,  4.0426
     &,  3.9867,  3.9661,  3.9345,  3.9200,  3.8883,  3.8498,  4.0496
     &,  3.8145,  4.0881,  3.8756,  3.7271,  4.3128,  4.5242,  4.3578
     &/)
      R0AB(1751:1820)=(/
     &   4.2870,  3.7796,  3.7318,  3.6364,  3.5854,  3.5726,  3.5378
     &,  3.5155,  4.0527,  4.0478,  4.2630,  4.0322,  4.3449,  4.1421
     &,  3.9975,  4.5499,  4.8825,  4.5601,  4.5950,  4.5702,  2.9046
     &,  3.2044,  3.5621,  3.8078,  3.6185,  3.6220,  3.4019,  3.4359
     &,  3.2110,  3.0635,  3.7037,  3.9958,  3.8792,  4.0194,  3.7460
     &,  3.9517,  3.7128,  3.5474,  4.0872,  4.3138,  4.1906,  4.1593
     &,  4.0973,  3.4919,  4.0216,  3.9657,  3.9454,  3.9134,  3.8986
     &,  3.8669,  3.8289,  4.0323,  3.7954,  4.0725,  3.8598,  3.7113
     &,  4.2896,  4.5021,  4.3325,  4.2645,  3.7571,  3.7083,  3.6136
     &,  3.5628,  3.5507,  3.5155,  3.4929,  4.0297,  4.0234,  4.2442
     &/)
      R0AB(1821:1890)=(/
     &   4.0112,  4.3274,  4.1240,  3.9793,  4.5257,  4.8568,  4.5353
     &,  4.5733,  4.5485,  4.5271,  2.8878,  3.1890,  3.5412,  3.7908
     &,  3.5974,  3.6078,  3.3871,  3.4243,  3.1992,  3.0513,  3.6831
     &,  3.9784,  3.8579,  4.0049,  3.7304,  3.9392,  3.7002,  3.5347
     &,  4.0657,  4.2955,  4.1705,  4.1424,  4.0800,  3.4717,  4.0043
     &,  3.9485,  3.9286,  3.8965,  3.8815,  3.8500,  3.8073,  4.0180
     &,  3.7796,  4.0598,  3.8470,  3.6983,  4.2678,  4.4830,  4.3132
     &,  4.2444,  3.7370,  3.6876,  3.5935,  3.5428,  3.5314,  3.4958
     &,  3.4730,  4.0117,  4.0043,  4.2287,  3.9939,  4.3134,  4.1096
     &,  3.9646,  4.5032,  4.8356,  4.5156,  4.5544,  4.5297,  4.5083
     &/)
      R0AB(1891:1960)=(/
     &   4.4896,  2.8709,  3.1737,  3.5199,  3.7734,  3.5802,  3.5934
     &,  3.3724,  3.4128,  3.1877,  3.0396,  3.6624,  3.9608,  3.8397
     &,  3.9893,  3.7145,  3.9266,  3.6877,  3.5222,  4.0448,  4.2771
     &,  4.1523,  4.1247,  4.0626,  3.4530,  3.9866,  3.9310,  3.9115
     &,  3.8792,  3.8641,  3.8326,  3.7892,  4.0025,  3.7636,  4.0471
     &,  3.8343,  3.6854,  4.2464,  4.4635,  4.2939,  4.2252,  3.7169
     &,  3.6675,  3.5739,  3.5235,  3.5126,  3.4768,  3.4537,  3.9932
     &,  3.9854,  4.2123,  3.9765,  4.2992,  4.0951,  3.9500,  4.4811
     &,  4.8135,  4.4959,  4.5351,  4.5105,  4.4891,  4.4705,  4.4515
     &,  2.8568,  3.1608,  3.5050,  3.7598,  3.5665,  3.5803,  3.3601
     &/)
      R0AB(1961:2030)=(/
     &   3.4031,  3.1779,  3.0296,  3.6479,  3.9471,  3.8262,  3.9773
     &,  3.7015,  3.9162,  3.6771,  3.5115,  4.0306,  4.2634,  4.1385
     &,  4.1116,  4.0489,  3.4366,  3.9732,  3.9176,  3.8983,  3.8659
     &,  3.8507,  3.8191,  3.7757,  3.9907,  3.7506,  4.0365,  3.8235
     &,  3.6745,  4.2314,  4.4490,  4.2792,  4.2105,  3.7003,  3.6510
     &,  3.5578,  3.5075,  3.4971,  3.4609,  3.4377,  3.9788,  3.9712
     &,  4.1997,  3.9624,  4.2877,  4.0831,  3.9378,  4.4655,  4.7974
     &,  4.4813,  4.5209,  4.4964,  4.4750,  4.4565,  4.4375,  4.4234
     &,  2.6798,  3.0151,  3.2586,  3.5292,  3.5391,  3.4902,  3.2887
     &,  3.3322,  3.1228,  2.9888,  3.4012,  3.7145,  3.7830,  3.6665
     &/)
      R0AB(2031:2100)=(/
     &   3.5898,  3.8077,  3.5810,  3.4265,  3.7726,  4.0307,  3.9763
     &,  3.8890,  3.8489,  3.2706,  3.7595,  3.6984,  3.6772,  3.6428
     &,  3.6243,  3.5951,  3.7497,  3.6775,  3.6364,  3.9203,  3.7157
     &,  3.5746,  3.9494,  4.2076,  4.1563,  4.0508,  3.5329,  3.4780
     &,  3.3731,  3.3126,  3.2846,  3.2426,  3.2135,  3.7491,  3.9006
     &,  3.8332,  3.8029,  4.1436,  3.9407,  3.7998,  4.1663,  4.5309
     &,  4.3481,  4.2911,  4.2671,  4.2415,  4.2230,  4.2047,  4.1908
     &,  4.1243,  2.5189,  2.9703,  3.3063,  3.6235,  3.4517,  3.3989
     &,  3.2107,  3.2434,  3.0094,  2.8580,  3.4253,  3.8157,  3.7258
     &,  3.6132,  3.5297,  3.7566,  3.5095,  3.3368,  3.7890,  4.1298
     &/)
      R0AB(2101:2170)=(/
     &   4.0190,  3.9573,  3.9237,  3.2677,  3.8480,  3.8157,  3.7656
     &,  3.7317,  3.7126,  3.6814,  3.6793,  3.6218,  3.5788,  3.8763
     &,  3.6572,  3.5022,  3.9737,  4.3255,  4.1828,  4.1158,  3.5078
     &,  3.4595,  3.3600,  3.3088,  3.2575,  3.2164,  3.1856,  3.8522
     &,  3.8665,  3.8075,  3.7772,  4.1391,  3.9296,  3.7772,  4.2134
     &,  4.7308,  4.3787,  4.3894,  4.3649,  4.3441,  4.3257,  4.3073
     &,  4.2941,  4.1252,  4.2427,  3.0481,  2.9584,  3.6919,  3.5990
     &,  3.8881,  3.4209,  3.1606,  3.1938,  2.9975,  2.8646,  3.8138
     &,  3.7935,  3.7081,  3.9155,  3.5910,  3.4808,  3.4886,  3.3397
     &,  4.1336,  4.1122,  3.9888,  3.9543,  3.8917,  3.5894,  3.8131
     &/)
      R0AB(2171:2240)=(/
     &   3.7635,  3.7419,  3.7071,  3.6880,  3.6574,  3.6546,  3.9375
     &,  3.6579,  3.5870,  3.6361,  3.5039,  4.3149,  4.2978,  4.1321
     &,  4.1298,  3.8164,  3.7680,  3.7154,  3.6858,  3.6709,  3.6666
     &,  3.6517,  3.8174,  3.8608,  4.1805,  3.9102,  3.8394,  3.8968
     &,  3.7673,  4.5274,  4.6682,  4.3344,  4.3639,  4.3384,  4.3162
     &,  4.2972,  4.2779,  4.2636,  4.0253,  4.1168,  4.1541,  2.8136
     &,  3.0951,  3.4635,  3.6875,  3.4987,  3.5183,  3.2937,  3.3580
     &,  3.1325,  2.9832,  3.6078,  3.8757,  3.7616,  3.9222,  3.6370
     &,  3.8647,  3.6256,  3.4595,  3.9874,  4.1938,  4.0679,  4.0430
     &,  3.9781,  3.3886,  3.9008,  3.8463,  3.8288,  3.7950,  3.7790
     &/)
      R0AB(2241:2310)=(/
     &   3.7472,  3.7117,  3.9371,  3.6873,  3.9846,  3.7709,  3.6210
     &,  4.1812,  4.3750,  4.2044,  4.1340,  3.6459,  3.5929,  3.5036
     &,  3.4577,  3.4528,  3.4146,  3.3904,  3.9014,  3.9031,  4.1443
     &,  3.8961,  4.2295,  4.0227,  3.8763,  4.4086,  4.7097,  4.4064
     &,  4.4488,  4.4243,  4.4029,  4.3842,  4.3655,  4.3514,  4.1162
     &,  4.2205,  4.1953,  4.2794,  2.8032,  3.0805,  3.4519,  3.6700
     &,  3.4827,  3.5050,  3.2799,  3.3482,  3.1233,  2.9747,  3.5971
     &,  3.8586,  3.7461,  3.9100,  3.6228,  3.8535,  3.6147,  3.4490
     &,  3.9764,  4.1773,  4.0511,  4.0270,  3.9614,  3.3754,  3.8836
     &,  3.8291,  3.8121,  3.7780,  3.7619,  3.7300,  3.6965,  3.9253
     &/)
      R0AB(2311:2380)=(/
     &   3.6734,  3.9733,  3.7597,  3.6099,  4.1683,  4.3572,  4.1862
     &,  4.1153,  3.6312,  3.5772,  3.4881,  3.4429,  3.4395,  3.4009
     &,  3.3766,  3.8827,  3.8868,  4.1316,  3.8807,  4.2164,  4.0092
     &,  3.8627,  4.3936,  4.6871,  4.3882,  4.4316,  4.4073,  4.3858
     &,  4.3672,  4.3485,  4.3344,  4.0984,  4.2036,  4.1791,  4.2622
     &,  4.2450,  2.7967,  3.0689,  3.4445,  3.6581,  3.4717,  3.4951
     &,  3.2694,  3.3397,  3.1147,  2.9661,  3.5898,  3.8468,  3.7358
     &,  3.9014,  3.6129,  3.8443,  3.6054,  3.4396,  3.9683,  4.1656
     &,  4.0394,  4.0158,  3.9498,  3.3677,  3.8718,  3.8164,  3.8005
     &,  3.7662,  3.7500,  3.7181,  3.6863,  3.9170,  3.6637,  3.9641
     &/)
      R0AB(2381:2450)=(/
     &   3.7503,  3.6004,  4.1590,  4.3448,  4.1739,  4.1029,  3.6224
     &,  3.5677,  3.4785,  3.4314,  3.4313,  3.3923,  3.3680,  3.8698
     &,  3.8758,  4.1229,  3.8704,  4.2063,  3.9987,  3.8519,  4.3832
     &,  4.6728,  4.3759,  4.4195,  4.3952,  4.3737,  4.3551,  4.3364
     &,  4.3223,  4.0861,  4.1911,  4.1676,  4.2501,  4.2329,  4.2208
     &,  2.7897,  3.0636,  3.4344,  3.6480,  3.4626,  3.4892,  3.2626
     &,  3.3344,  3.1088,  2.9597,  3.5804,  3.8359,  3.7251,  3.8940
     &,  3.6047,  3.8375,  3.5990,  3.4329,  3.9597,  4.1542,  4.0278
     &,  4.0048,  3.9390,  3.3571,  3.8608,  3.8056,  3.7899,  3.7560
     &,  3.7400,  3.7081,  3.6758,  3.9095,  3.6552,  3.9572,  3.7436
     &/)
      R0AB(2451:2520)=(/
     &   3.5933,  4.1508,  4.3337,  4.1624,  4.0916,  3.6126,  3.5582
     &,  3.4684,  3.4212,  3.4207,  3.3829,  3.3586,  3.8604,  3.8658
     &,  4.1156,  3.8620,  4.1994,  3.9917,  3.8446,  4.3750,  4.6617
     &,  4.3644,  4.4083,  4.3840,  4.3625,  4.3439,  4.3253,  4.3112
     &,  4.0745,  4.1807,  4.1578,  4.2390,  4.2218,  4.2097,  4.1986
     &,  2.8395,  3.0081,  3.3171,  3.4878,  3.5360,  3.5145,  3.2809
     &,  3.3307,  3.1260,  2.9940,  3.4741,  3.6675,  3.7832,  3.6787
     &,  3.6156,  3.8041,  3.5813,  3.4301,  3.8480,  3.9849,  3.9314
     &,  3.8405,  3.8029,  3.2962,  3.7104,  3.6515,  3.6378,  3.6020
     &,  3.5849,  3.5550,  3.7494,  3.6893,  3.6666,  3.9170,  3.7150
     &/)
      R0AB(2521:2590)=(/
     &   3.5760,  4.0268,  4.1596,  4.1107,  3.9995,  3.5574,  3.5103
     &,  3.4163,  3.3655,  3.3677,  3.3243,  3.2975,  3.7071,  3.9047
     &,  3.8514,  3.8422,  3.8022,  3.9323,  3.7932,  4.2343,  4.4583
     &,  4.3115,  4.2457,  4.2213,  4.1945,  4.1756,  4.1569,  4.1424
     &,  4.0620,  4.0494,  3.9953,  4.0694,  4.0516,  4.0396,  4.0280
     &,  4.0130,  2.9007,  2.9674,  3.8174,  3.5856,  3.6486,  3.5339
     &,  3.2832,  3.3154,  3.1144,  2.9866,  3.9618,  3.8430,  3.9980
     &,  3.8134,  3.6652,  3.7985,  3.5756,  3.4207,  4.4061,  4.2817
     &,  4.1477,  4.0616,  3.9979,  3.6492,  3.8833,  3.8027,  3.7660
     &,  3.7183,  3.6954,  3.6525,  3.9669,  3.8371,  3.7325,  3.9160
     &/)
      R0AB(2591:2660)=(/
     &   3.7156,  3.5714,  4.6036,  4.4620,  4.3092,  4.2122,  3.8478
     &,  3.7572,  3.6597,  3.5969,  3.5575,  3.5386,  3.5153,  3.7818
     &,  4.1335,  4.0153,  3.9177,  3.8603,  3.9365,  3.7906,  4.7936
     &,  4.7410,  4.5461,  4.5662,  4.5340,  4.5059,  4.4832,  4.4604
     &,  4.4429,  4.2346,  4.4204,  4.3119,  4.3450,  4.3193,  4.3035
     &,  4.2933,  4.1582,  4.2450,  2.8559,  2.9050,  3.8325,  3.5442
     &,  3.5077,  3.4905,  3.2396,  3.2720,  3.0726,  2.9467,  3.9644
     &,  3.8050,  3.8981,  3.7762,  3.6216,  3.7531,  3.5297,  3.3742
     &,  4.3814,  4.2818,  4.1026,  4.0294,  3.9640,  3.6208,  3.8464
     &,  3.7648,  3.7281,  3.6790,  3.6542,  3.6117,  3.8650,  3.8010
     &/)
      R0AB(2661:2730)=(/
     &   3.6894,  3.8713,  3.6699,  3.5244,  4.5151,  4.4517,  4.2538
     &,  4.1483,  3.8641,  3.7244,  3.6243,  3.5589,  3.5172,  3.4973
     &,  3.4715,  3.7340,  4.0316,  3.9958,  3.8687,  3.8115,  3.8862
     &,  3.7379,  4.7091,  4.7156,  4.5199,  4.5542,  4.5230,  4.4959
     &,  4.4750,  4.4529,  4.4361,  4.1774,  4.3774,  4.2963,  4.3406
     &,  4.3159,  4.3006,  4.2910,  4.1008,  4.1568,  4.0980,  2.8110
     &,  2.8520,  3.7480,  3.5105,  3.4346,  3.3461,  3.1971,  3.2326
     &,  3.0329,  2.9070,  3.8823,  3.7928,  3.8264,  3.7006,  3.5797
     &,  3.7141,  3.4894,  3.3326,  4.3048,  4.2217,  4.0786,  3.9900
     &,  3.9357,  3.6331,  3.8333,  3.7317,  3.6957,  3.6460,  3.6197
     &/)
      R0AB(2731:2800)=(/
     &   3.5779,  3.7909,  3.7257,  3.6476,  3.5729,  3.6304,  3.4834
     &,  4.4368,  4.3921,  4.2207,  4.1133,  3.8067,  3.7421,  3.6140
     &,  3.5491,  3.5077,  3.4887,  3.4623,  3.6956,  3.9568,  3.8976
     &,  3.8240,  3.7684,  3.8451,  3.6949,  4.6318,  4.6559,  4.4533
     &,  4.4956,  4.4641,  4.4366,  4.4155,  4.3936,  4.3764,  4.1302
     &,  4.3398,  4.2283,  4.2796,  4.2547,  4.2391,  4.2296,  4.0699
     &,  4.1083,  4.0319,  3.9855,  2.7676,  2.8078,  3.6725,  3.4804
     &,  3.3775,  3.2411,  3.1581,  3.1983,  2.9973,  2.8705,  3.8070
     &,  3.7392,  3.7668,  3.6263,  3.5402,  3.6807,  3.4545,  3.2962
     &,  4.2283,  4.1698,  4.0240,  3.9341,  3.8711,  3.5489,  3.7798
     &/)
      R0AB(2801:2870)=(/
     &   3.7000,  3.6654,  3.6154,  3.5882,  3.5472,  3.7289,  3.6510
     &,  3.6078,  3.5355,  3.5963,  3.4480,  4.3587,  4.3390,  4.1635
     &,  4.0536,  3.7193,  3.6529,  3.5512,  3.4837,  3.4400,  3.4191
     &,  3.3891,  3.6622,  3.8934,  3.8235,  3.7823,  3.7292,  3.8106
     &,  3.6589,  4.5535,  4.6013,  4.3961,  4.4423,  4.4109,  4.3835
     &,  4.3625,  4.3407,  4.3237,  4.0863,  4.2835,  4.1675,  4.2272
     &,  4.2025,  4.1869,  4.1774,  4.0126,  4.0460,  3.9815,  3.9340
     &,  3.8955,  2.6912,  2.7604,  3.6037,  3.4194,  3.3094,  3.1710
     &,  3.0862,  3.1789,  2.9738,  2.8427,  3.7378,  3.6742,  3.6928
     &,  3.5512,  3.4614,  3.4087,  3.4201,  3.2607,  4.1527,  4.0977
     &/)
      R0AB(2871:2940)=(/
     &   3.9523,  3.8628,  3.8002,  3.4759,  3.7102,  3.6466,  3.6106
     &,  3.5580,  3.5282,  3.4878,  3.6547,  3.5763,  3.5289,  3.5086
     &,  3.5593,  3.4099,  4.2788,  4.2624,  4.0873,  3.9770,  3.6407
     &,  3.5743,  3.5178,  3.4753,  3.3931,  3.3694,  3.3339,  3.6002
     &,  3.8164,  3.7478,  3.7028,  3.6952,  3.7669,  3.6137,  4.4698
     &,  4.5488,  4.3168,  4.3646,  4.3338,  4.3067,  4.2860,  4.2645
     &,  4.2478,  4.0067,  4.2349,  4.0958,  4.1543,  4.1302,  4.1141
     &,  4.1048,  3.9410,  3.9595,  3.8941,  3.8465,  3.8089,  3.7490
     &,  2.7895,  2.5849,  3.6484,  3.0162,  3.1267,  3.2125,  3.0043
     &,  2.9572,  2.8197,  2.7261,  3.7701,  3.2446,  3.5239,  3.4696
     &/)
      R0AB(2941:3010)=(/
     &   3.4261,  3.3508,  3.1968,  3.0848,  4.1496,  3.6598,  3.5111
     &,  3.4199,  3.3809,  3.5382,  3.2572,  3.2100,  3.1917,  3.1519
     &,  3.1198,  3.1005,  3.5071,  3.5086,  3.5073,  3.4509,  3.3120
     &,  3.2082,  4.2611,  3.8117,  3.6988,  3.5646,  3.6925,  3.6295
     &,  3.5383,  3.4910,  3.4625,  3.4233,  3.4007,  3.2329,  3.6723
     &,  3.6845,  3.6876,  3.6197,  3.4799,  3.3737,  4.4341,  4.0525
     &,  3.9011,  3.8945,  3.8635,  3.8368,  3.8153,  3.7936,  3.7758
     &,  3.4944,  3.4873,  3.9040,  3.7110,  3.6922,  3.6799,  3.6724
     &,  3.5622,  3.6081,  3.5426,  3.4922,  3.4498,  3.3984,  3.4456
     &,  2.7522,  2.5524,  3.5742,  2.9508,  3.0751,  3.0158,  2.9644
     &/)
      R0AB(3011:3080)=(/
     &   2.8338,  2.7891,  2.6933,  3.6926,  3.1814,  3.4528,  3.4186
     &,  3.3836,  3.2213,  3.1626,  3.0507,  4.0548,  3.5312,  3.4244
     &,  3.3409,  3.2810,  3.4782,  3.1905,  3.1494,  3.1221,  3.1128
     &,  3.0853,  3.0384,  3.4366,  3.4562,  3.4638,  3.3211,  3.2762
     &,  3.1730,  4.1632,  3.6825,  3.5822,  3.4870,  3.6325,  3.5740
     &,  3.4733,  3.4247,  3.3969,  3.3764,  3.3525,  3.1984,  3.5989
     &,  3.6299,  3.6433,  3.4937,  3.4417,  3.3365,  4.3304,  3.9242
     &,  3.7793,  3.7623,  3.7327,  3.7071,  3.6860,  3.6650,  3.6476
     &,  3.3849,  3.3534,  3.8216,  3.5870,  3.5695,  3.5584,  3.5508
     &,  3.4856,  3.5523,  3.4934,  3.4464,  3.4055,  3.3551,  3.3888
     &/)
      R0AB(3081:3150)=(/
     &   3.3525,  2.7202,  2.5183,  3.4947,  2.8731,  3.0198,  3.1457
     &,  2.9276,  2.7826,  2.7574,  2.6606,  3.6090,  3.0581,  3.3747
     &,  3.3677,  3.3450,  3.1651,  3.1259,  3.0147,  3.9498,  3.3857
     &,  3.2917,  3.2154,  3.1604,  3.4174,  3.0735,  3.0342,  3.0096
     &,  3.0136,  2.9855,  2.9680,  3.3604,  3.4037,  3.4243,  3.2633
     &,  3.1810,  3.1351,  4.0557,  3.5368,  3.4526,  3.3699,  3.5707
     &,  3.5184,  3.4085,  3.3595,  3.3333,  3.3143,  3.3041,  3.1094
     &,  3.5193,  3.5745,  3.6025,  3.4338,  3.3448,  3.2952,  4.2158
     &,  3.7802,  3.6431,  3.6129,  3.5853,  3.5610,  3.5406,  3.5204
     &,  3.5036,  3.2679,  3.2162,  3.7068,  3.4483,  3.4323,  3.4221
     &/)
      R0AB(3151:3220)=(/
     &   3.4138,  3.3652,  3.4576,  3.4053,  3.3618,  3.3224,  3.2711
     &,  3.3326,  3.2950,  3.2564,  2.5315,  2.6104,  3.2734,  3.2299
     &,  3.1090,  2.9942,  2.9159,  2.8324,  2.8350,  2.7216,  3.3994
     &,  3.4475,  3.4354,  3.3438,  3.2807,  3.2169,  3.2677,  3.1296
     &,  3.7493,  3.8075,  3.6846,  3.6104,  3.5577,  3.2052,  3.4803
     &,  3.4236,  3.3845,  3.3640,  3.3365,  3.3010,  3.3938,  3.3624
     &,  3.3440,  3.3132,  3.4035,  3.2754,  3.8701,  3.9523,  3.8018
     &,  3.7149,  3.3673,  3.3199,  3.2483,  3.2069,  3.1793,  3.1558
     &,  3.1395,  3.4097,  3.5410,  3.5228,  3.5116,  3.4921,  3.4781
     &,  3.4690,  4.0420,  4.1759,  4.0078,  4.0450,  4.0189,  3.9952
     &/)
      R0AB(3221:3290)=(/
     &   3.9770,  3.9583,  3.9434,  3.7217,  3.8228,  3.7826,  3.8640
     &,  3.8446,  3.8314,  3.8225,  3.6817,  3.7068,  3.6555,  3.6159
     &,  3.5831,  3.5257,  3.2133,  3.1689,  3.1196,  3.3599,  2.9852
     &,  2.7881,  3.5284,  3.3493,  3.6958,  3.3642,  3.1568,  3.0055
     &,  2.9558,  2.8393,  3.6287,  3.5283,  4.1511,  3.8259,  3.6066
     &,  3.4527,  3.3480,  3.2713,  3.9037,  3.8361,  3.8579,  3.7311
     &,  3.6575,  3.5176,  3.5693,  3.5157,  3.4814,  3.4559,  3.4445
     &,  3.4160,  4.1231,  3.8543,  3.6816,  3.5602,  3.4798,  3.4208
     &,  4.0542,  4.0139,  4.0165,  3.9412,  3.7698,  3.6915,  3.6043
     &,  3.5639,  3.5416,  3.5247,  3.5153,  3.5654,  4.2862,  4.0437
     &/)
      R0AB(3291:3360)=(/
     &   3.8871,  3.7741,  3.6985,  3.6413,  4.2345,  4.3663,  4.3257
     &,  4.0869,  4.0612,  4.0364,  4.0170,  3.9978,  3.9834,  3.9137
     &,  3.8825,  3.8758,  3.9143,  3.8976,  3.8864,  3.8768,  3.9190
     &,  4.1613,  4.0566,  3.9784,  3.9116,  3.8326,  3.7122,  3.6378
     &,  3.5576,  3.5457,  4.3127,  3.1160,  2.8482,  4.0739,  3.3599
     &,  3.5698,  3.5366,  3.2854,  3.1039,  2.9953,  2.9192,  4.1432
     &,  3.5320,  3.9478,  4.0231,  3.7509,  3.5604,  3.4340,  3.3426
     &,  4.3328,  3.8288,  3.7822,  3.7909,  3.6907,  3.6864,  3.5793
     &,  3.5221,  3.4883,  3.4649,  3.4514,  3.4301,  3.9256,  4.0596
     &,  3.8307,  3.6702,  3.5651,  3.4884,  4.4182,  4.2516,  3.9687
     &/)
      R0AB(3361:3430)=(/
     &   3.9186,  3.9485,  3.8370,  3.7255,  3.6744,  3.6476,  3.6295
     &,  3.6193,  3.5659,  4.0663,  4.2309,  4.0183,  3.8680,  3.7672
     &,  3.6923,  4.5240,  4.4834,  4.1570,  4.3204,  4.2993,  4.2804
     &,  4.2647,  4.2481,  4.2354,  3.8626,  3.8448,  4.2267,  4.1799
     &,  4.1670,  3.8738,  3.8643,  3.8796,  4.0575,  4.0354,  3.9365
     &,  3.8611,  3.7847,  3.7388,  3.6826,  3.6251,  3.5492,  4.0889
     &,  4.2764,  3.1416,  2.8325,  3.7735,  3.3787,  3.4632,  3.5923
     &,  3.3214,  3.1285,  3.0147,  2.9366,  3.8527,  3.5602,  3.8131
     &,  3.8349,  3.7995,  3.5919,  3.4539,  3.3540,  4.0654,  3.8603
     &,  3.7972,  3.7358,  3.7392,  3.8157,  3.6055,  3.5438,  3.5089
     &/)
      R0AB(3431:3500)=(/
     &   3.4853,  3.4698,  3.4508,  3.7882,  3.8682,  3.8837,  3.7055
     &,  3.5870,  3.5000,  4.1573,  4.0005,  3.9568,  3.8936,  3.9990
     &,  3.9433,  3.8172,  3.7566,  3.7246,  3.7033,  3.6900,  3.5697
     &,  3.9183,  4.0262,  4.0659,  3.8969,  3.7809,  3.6949,  4.2765
     &,  4.2312,  4.1401,  4.0815,  4.0580,  4.0369,  4.0194,  4.0017
     &,  3.9874,  3.8312,  3.8120,  3.9454,  3.9210,  3.9055,  3.8951
     &,  3.8866,  3.8689,  3.9603,  3.9109,  3.9122,  3.8233,  3.7438
     &,  3.7436,  3.6981,  3.6555,  3.5452,  3.9327,  4.0658,  4.1175
     &,  2.9664,  2.8209,  3.5547,  3.3796,  3.3985,  3.3164,  3.2364
     &,  3.1956,  3.0370,  2.9313,  3.6425,  3.5565,  3.7209,  3.7108
     &/)
      R0AB(3501:3570)=(/
     &   3.6639,  3.6484,  3.4745,  3.3492,  3.8755,  4.2457,  3.7758
     &,  3.7161,  3.6693,  3.6155,  3.5941,  3.5643,  3.5292,  3.4950
     &,  3.4720,  3.4503,  3.6936,  3.7392,  3.7388,  3.7602,  3.6078
     &,  3.4960,  3.9800,  4.3518,  4.2802,  3.8580,  3.8056,  3.7527
     &,  3.7019,  3.6615,  3.5768,  3.5330,  3.5038,  3.5639,  3.8192
     &,  3.8883,  3.9092,  3.9478,  3.7995,  3.6896,  4.1165,  4.5232
     &,  4.4357,  4.4226,  4.4031,  4.3860,  4.3721,  4.3580,  4.3466
     &,  4.2036,  4.2037,  3.8867,  4.2895,  4.2766,  4.2662,  4.2598
     &,  3.8408,  3.9169,  3.8681,  3.8250,  3.7855,  3.7501,  3.6753
     &,  3.5499,  3.4872,  3.5401,  3.8288,  3.9217,  3.9538,  4.0054
     &/)
      R0AB(3571:3640)=(/
     &   2.8388,  2.7890,  3.4329,  3.5593,  3.3488,  3.2486,  3.1615
     &,  3.1000,  3.0394,  2.9165,  3.5267,  3.7479,  3.6650,  3.6263
     &,  3.5658,  3.5224,  3.4762,  3.3342,  3.7738,  4.0333,  3.9568
     &,  3.8975,  3.8521,  3.4929,  3.7830,  3.7409,  3.7062,  3.6786
     &,  3.6471,  3.6208,  3.6337,  3.6519,  3.6363,  3.6278,  3.6110
     &,  3.4825,  3.8795,  4.1448,  4.0736,  4.0045,  3.6843,  3.6291
     &,  3.5741,  3.5312,  3.4974,  3.4472,  3.4034,  3.7131,  3.7557
     &,  3.7966,  3.8005,  3.8068,  3.8015,  3.6747,  4.0222,  4.3207
     &,  4.2347,  4.2191,  4.1990,  4.1811,  4.1666,  4.1521,  4.1401
     &,  3.9970,  3.9943,  3.9592,  4.0800,  4.0664,  4.0559,  4.0488
     &/)
      R0AB(3641:3710)=(/
     &   3.9882,  4.0035,  3.9539,  3.9138,  3.8798,  3.8355,  3.5359
     &,  3.4954,  3.3962,  3.5339,  3.7595,  3.8250,  3.8408,  3.8600
     &,  3.8644,  2.7412,  2.7489,  3.3374,  3.3950,  3.3076,  3.1910
     &,  3.0961,  3.0175,  3.0280,  2.8929,  3.4328,  3.5883,  3.6227
     &,  3.5616,  3.4894,  3.4241,  3.3641,  3.3120,  3.6815,  3.8789
     &,  3.8031,  3.7413,  3.6939,  3.4010,  3.6225,  3.5797,  3.5443
     &,  3.5139,  3.4923,  3.4642,  3.5860,  3.5849,  3.5570,  3.5257
     &,  3.4936,  3.4628,  3.7874,  3.9916,  3.9249,  3.8530,  3.5932
     &,  3.5355,  3.4757,  3.4306,  3.3953,  3.3646,  3.3390,  3.5637
     &,  3.7053,  3.7266,  3.7177,  3.6996,  3.6775,  3.6558,  3.9331
     &/)
      R0AB(3711:3780)=(/
     &   4.1655,  4.0879,  4.0681,  4.0479,  4.0299,  4.0152,  4.0006
     &,  3.9883,  3.8500,  3.8359,  3.8249,  3.9269,  3.9133,  3.9025
     &,  3.8948,  3.8422,  3.8509,  3.7990,  3.7570,  3.7219,  3.6762
     &,  3.4260,  3.3866,  3.3425,  3.5294,  3.7022,  3.7497,  3.7542
     &,  3.7494,  3.7370,  3.7216,  3.4155,  3.0522,  4.2541,  3.8218
     &,  4.0438,  3.5875,  3.3286,  3.1682,  3.0566,  2.9746,  4.3627
     &,  4.0249,  4.6947,  4.1718,  3.8639,  3.6735,  3.5435,  3.4479
     &,  4.6806,  4.3485,  4.2668,  4.1690,  4.1061,  4.1245,  4.0206
     &,  3.9765,  3.9458,  3.9217,  3.9075,  3.8813,  3.9947,  4.1989
     &,  3.9507,  3.7960,  3.6925,  3.6150,  4.8535,  4.5642,  4.4134
     &/)
      R0AB(3781:3850)=(/
     &   4.3688,  4.3396,  4.2879,  4.2166,  4.1888,  4.1768,  4.1660
     &,  4.1608,  4.0745,  4.2289,  4.4863,  4.2513,  4.0897,  3.9876
     &,  3.9061,  5.0690,  5.0446,  4.6186,  4.6078,  4.5780,  4.5538
     &,  4.5319,  4.5101,  4.4945,  4.1912,  4.2315,  4.5534,  4.4373
     &,  4.4224,  4.4120,  4.4040,  4.2634,  4.7770,  4.6890,  4.6107
     &,  4.5331,  4.4496,  4.4082,  4.3095,  4.2023,  4.0501,  4.2595
     &,  4.5497,  4.3056,  4.1506,  4.0574,  3.9725,  5.0796,  3.0548
     &,  3.3206,  3.8132,  3.9720,  3.7675,  3.7351,  3.5167,  3.5274
     &,  3.3085,  3.1653,  3.9500,  4.1730,  4.0613,  4.1493,  3.8823
     &,  4.0537,  3.8200,  3.6582,  4.3422,  4.5111,  4.3795,  4.3362
     &/)
      R0AB(3851:3920)=(/
     &   4.2751,  3.7103,  4.1973,  4.1385,  4.1129,  4.0800,  4.0647
     &,  4.0308,  4.0096,  4.1619,  3.9360,  4.1766,  3.9705,  3.8262
     &,  4.5348,  4.7025,  4.5268,  4.5076,  3.9562,  3.9065,  3.8119
     &,  3.7605,  3.7447,  3.7119,  3.6916,  4.1950,  4.2110,  4.3843
     &,  4.1631,  4.4427,  4.2463,  4.1054,  4.7693,  5.0649,  4.7365
     &,  4.7761,  4.7498,  4.7272,  4.7076,  4.6877,  4.6730,  4.4274
     &,  4.5473,  4.5169,  4.5975,  4.5793,  4.5667,  4.5559,  4.3804
     &,  4.6920,  4.6731,  4.6142,  4.5600,  4.4801,  4.0149,  3.8856
     &,  3.7407,  4.1545,  4.2253,  4.4229,  4.1923,  4.5022,  4.3059
     &,  4.1591,  4.7883,  4.9294,  3.3850,  3.4208,  3.7004,  3.8800
     &/)
      R0AB(3921:3990)=(/
     &   3.9886,  3.9040,  3.6719,  3.6547,  3.4625,  3.3370,  3.8394
     &,  4.0335,  4.2373,  4.3023,  4.0306,  4.1408,  3.9297,  3.7857
     &,  4.1907,  4.3230,  4.2664,  4.2173,  4.1482,  3.6823,  4.0711
     &,  4.0180,  4.0017,  3.9747,  3.9634,  3.9383,  4.1993,  4.3205
     &,  4.0821,  4.2547,  4.0659,  3.9359,  4.3952,  4.5176,  4.3888
     &,  4.3607,  3.9583,  3.9280,  3.8390,  3.7971,  3.7955,  3.7674
     &,  3.7521,  4.1062,  4.3633,  4.2991,  4.2767,  4.4857,  4.3039
     &,  4.1762,  4.6197,  4.8654,  4.6633,  4.5878,  4.5640,  4.5422
     &,  4.5231,  4.5042,  4.4901,  4.3282,  4.3978,  4.3483,  4.4202
     &,  4.4039,  4.3926,  4.3807,  4.2649,  4.6135,  4.5605,  4.5232
     &/)
      R0AB(3991:4060)=(/
     &   4.4676,  4.3948,  4.0989,  3.9864,  3.8596,  4.0942,  4.2720
     &,  4.3270,  4.3022,  4.5410,  4.3576,  4.2235,  4.6545,  4.7447
     &,  4.7043,  3.0942,  3.2075,  3.5152,  3.6659,  3.8289,  3.7459
     &,  3.5156,  3.5197,  3.3290,  3.2069,  3.6702,  3.8448,  4.0340
     &,  3.9509,  3.8585,  3.9894,  3.7787,  3.6365,  4.1425,  4.1618
     &,  4.0940,  4.0466,  3.9941,  3.5426,  3.8952,  3.8327,  3.8126
     &,  3.7796,  3.7635,  3.7356,  4.0047,  3.9655,  3.9116,  4.1010
     &,  3.9102,  3.7800,  4.2964,  4.3330,  4.2622,  4.2254,  3.8195
     &,  3.7560,  3.6513,  3.5941,  3.5810,  3.5420,  3.5178,  3.8861
     &,  4.1459,  4.1147,  4.0772,  4.3120,  4.1207,  3.9900,  4.4733
     &/)
      R0AB(4061:4130)=(/
     &   4.6157,  4.4580,  4.4194,  4.3954,  4.3739,  4.3531,  4.3343
     &,  4.3196,  4.2140,  4.2339,  4.1738,  4.2458,  4.2278,  4.2158
     &,  4.2039,  4.1658,  4.3595,  4.2857,  4.2444,  4.1855,  4.1122
     &,  3.7839,  3.6879,  3.5816,  3.8633,  4.1585,  4.1402,  4.1036
     &,  4.3694,  4.1735,  4.0368,  4.5095,  4.5538,  4.5240,  4.4252
     &,  3.0187,  3.1918,  3.5127,  3.6875,  3.7404,  3.6943,  3.4702
     &,  3.4888,  3.2914,  3.1643,  3.6669,  3.8724,  3.9940,  4.0816
     &,  3.8054,  3.9661,  3.7492,  3.6024,  4.0428,  4.1951,  4.1466
     &,  4.0515,  4.0075,  3.5020,  3.9158,  3.8546,  3.8342,  3.8008
     &,  3.7845,  3.7549,  3.9602,  3.8872,  3.8564,  4.0793,  3.8835
     &/)
      R0AB(4131:4200)=(/
     &   3.7495,  4.2213,  4.3704,  4.3300,  4.2121,  3.7643,  3.7130
     &,  3.6144,  3.5599,  3.5474,  3.5093,  3.4853,  3.9075,  4.1115
     &,  4.0473,  4.0318,  4.2999,  4.1050,  3.9710,  4.4320,  4.6706
     &,  4.5273,  4.4581,  4.4332,  4.4064,  4.3873,  4.3684,  4.3537
     &,  4.2728,  4.2549,  4.2032,  4.2794,  4.2613,  4.2491,  4.2375
     &,  4.2322,  4.3665,  4.3061,  4.2714,  4.2155,  4.1416,  3.7660
     &,  3.6628,  3.5476,  3.8790,  4.1233,  4.0738,  4.0575,  4.3575
     &,  4.1586,  4.0183,  4.4593,  4.5927,  4.4865,  4.3813,  4.4594
     &,  2.9875,  3.1674,  3.4971,  3.6715,  3.7114,  3.6692,  3.4446
     &,  3.4676,  3.2685,  3.1405,  3.6546,  3.8579,  3.9637,  4.0581
     &/)
      R0AB(4201:4270)=(/
     &   3.7796,  3.9463,  3.7275,  3.5792,  4.0295,  4.1824,  4.1247
     &,  4.0357,  3.9926,  3.4827,  3.9007,  3.8392,  3.8191,  3.7851
     &,  3.7687,  3.7387,  3.9290,  3.8606,  3.8306,  4.0601,  3.8625
     &,  3.7269,  4.2062,  4.3566,  4.3022,  4.1929,  3.7401,  3.6888
     &,  3.5900,  3.5350,  3.5226,  3.4838,  3.4594,  3.8888,  4.0813
     &,  4.0209,  4.0059,  4.2810,  4.0843,  3.9486,  4.4162,  4.6542
     &,  4.5005,  4.4444,  4.4196,  4.3933,  4.3741,  4.3552,  4.3406
     &,  4.2484,  4.2413,  4.1907,  4.2656,  4.2474,  4.2352,  4.2236
     &,  4.2068,  4.3410,  4.2817,  4.2479,  4.1921,  4.1182,  3.7346
     &,  3.6314,  3.5168,  3.8582,  4.0927,  4.0469,  4.0313,  4.3391
     &/)
      R0AB(4271:4340)=(/
     &   4.1381,  3.9962,  4.4429,  4.5787,  4.4731,  4.3588,  4.4270
     &,  4.3957,  2.9659,  3.1442,  3.4795,  3.6503,  3.6814,  3.6476
     &,  3.4222,  3.4491,  3.2494,  3.1209,  3.6324,  3.8375,  3.9397
     &,  3.8311,  3.7581,  3.9274,  3.7085,  3.5598,  4.0080,  4.1641
     &,  4.1057,  4.0158,  3.9726,  3.4667,  3.8802,  3.8188,  3.7989
     &,  3.7644,  3.7474,  3.7173,  3.9049,  3.8424,  3.8095,  4.0412
     &,  3.8436,  3.7077,  4.1837,  4.3366,  4.2816,  4.1686,  3.7293
     &,  3.6709,  3.5700,  3.5153,  3.5039,  3.4684,  3.4437,  3.8663
     &,  4.0575,  4.0020,  3.9842,  4.2612,  4.0643,  3.9285,  4.3928
     &,  4.6308,  4.4799,  4.4244,  4.3996,  4.3737,  4.3547,  4.3358
     &/)
      R0AB(4341:4410)=(/
     &   4.3212,  4.2275,  4.2216,  4.1676,  4.2465,  4.2283,  4.2161
     &,  4.2045,  4.1841,  4.3135,  4.2562,  4.2226,  4.1667,  4.0932
     &,  3.7134,  3.6109,  3.4962,  3.8352,  4.0688,  4.0281,  4.0099
     &,  4.3199,  4.1188,  3.9768,  4.4192,  4.5577,  4.4516,  4.3365
     &,  4.4058,  4.3745,  4.3539,  2.8763,  3.1294,  3.5598,  3.7465
     &,  3.5659,  3.5816,  3.3599,  3.4024,  3.1877,  3.0484,  3.7009
     &,  3.9451,  3.8465,  3.9873,  3.7079,  3.9083,  3.6756,  3.5150
     &,  4.0829,  4.2780,  4.1511,  4.1260,  4.0571,  3.4865,  3.9744
     &,  3.9150,  3.8930,  3.8578,  3.8402,  3.8073,  3.7977,  4.0036
     &,  3.7604,  4.0288,  3.8210,  3.6757,  4.2646,  4.4558,  4.2862
     &/)
      R0AB(4411:4465)=(/
     &   4.2122,  3.7088,  3.6729,  3.5800,  3.5276,  3.5165,  3.4783
     &,  3.4539,  3.9553,  3.9818,  4.2040,  3.9604,  4.2718,  4.0689
     &,  3.9253,  4.4869,  4.7792,  4.4918,  4.5342,  4.5090,  4.4868
     &,  4.4680,  4.4486,  4.4341,  4.2023,  4.3122,  4.2710,  4.3587
     &,  4.3407,  4.3281,  4.3174,  4.1499,  4.3940,  4.3895,  4.3260
     &,  4.2725,  4.1961,  3.7361,  3.6193,  3.4916,  3.9115,  3.9914
     &,  3.9809,  3.9866,  4.3329,  4.1276,  3.9782,  4.5097,  4.6769
     &,  4.5158,  4.3291,  4.3609,  4.3462,  4.3265,  4.4341
     &/)

      K=0
      DO I=1,N_ELEM
         DO J=1,I
            K=K+1
            R(I,J)=R0AB(K)*ANG_TO_AU
            R(J,I)=R0AB(K)*ANG_TO_AU
         ENDDO
      ENDDO

      RETURN
      END



      SUBROUTINE SET_R0(N_ELEM_D2,R00)
      implicit none
      !Pass in as R00 space for R0
      !Load from data statement 
      DOUBLE PRECISION R0(86)
      INTEGER N_ELEM_D2,I
      DOUBLE PRECISION R00(N_ELEM_D2)
c the published radii in S.Grimme, J.Comput.Chem. 27, (2006), 1787-1799 (tab 1)
c refer to the following values multiplied by 1.1 (rs6 in this code)
C                        1      2
      DATA      R0/ 0.91d0,0.92d0,
C             3       4     5      6      7      8       9    10
     &      0.75d0,1.28d0,1.35d0,1.32d0,1.27d0,1.22d0,1.17d0,1.13d0,
C             11      12    13     14     15     16      17   18
     &      1.04d0,1.24d0,1.49d0,1.56d0,1.55d0,1.53d0,1.49d0,1.45d0,
C             19      20
     &      1.35d0,1.34d0,
C             21      22    23     24     25
     &      1.42d0,1.42d0,1.42d0,1.42d0,1.42d0,
C             26      27    28     29     30
     &      1.42d0,1.42d0,1.42d0,1.42d0,1.42d0,
C             31      32    33     34     35     36
     &      1.50d0,1.57d0,1.60d0,1.61d0,1.59d0,1.57d0,
C             37      38
     &      1.48d0,1.46d0,
C             39      40    41     42     43
     &      1.49d0,1.49d0,1.49d0,1.49d0,1.49d0,
C             44      45    46     47     48
     &      1.49d0,1.49d0,1.49d0,1.49d0,1.49d0,
C             49      50    51     52     53     54
     &      1.52d0,1.64d0,1.71d0,1.72d0,1.72d0,1.71d0,
C             55      56      57      58      59      60      61   
     &      1.638d0,1.602d0,1.564d0,1.594d0,1.594d0,1.594d0,1.594d0,
C             62      63      64      65      66      67      68   
     &      1.594d0,1.594d0,1.594d0,1.594d0,1.594d0,1.594d0,1.594d0,
C             69      70      71
     &      1.594d0,1.594d0,1.594d0,
C             72      73      74      75      76      77      78   
     &      1.625d0,1.611d0,1.611d0,1.611d0,1.611d0,1.611d0,1.611d0,
C             79
     &      1.611d0,
C             80      81      82      83      84      85      86   
     &      1.598d0,1.805d0,1.767d0,1.725d0,1.823d0,1.810d0,1.749d0/
      
      DO I = 1,N_ELEM_D2
        R00(I) = R0(I)
      ENDDO

      RETURN
      END


      SUBROUTINE SET_C6(N_ELEM_D2,C600)
      implicit none
      !Pass in as C600 space for C6
      !Load from data statement 
      DOUBLE PRECISION C6(86)
      INTEGER N_ELEM_D2,I
      DOUBLE PRECISION C600(N_ELEM_D2)
C                      1      2
      DATA      C6/0.14d0,0.08d0,
C          3       4     5      6      7      8       9    10
     &   1.61d0,1.61d0,3.13d0,1.75d0,1.23d0,0.70d0,0.75d0,0.63d0,
C          11      12    13     14     15     16      17   18
     &   5.71d0,5.71d0,10.79d0,9.23d0,7.84d0,5.57d0,5.07d0,4.61d0,
C          19      20    21     22     23
     &   10.8d0,10.8d0,10.8d0,10.8d0,10.8d0,
C          24      25    26     27     28     29     30     31
     &   10.8d0,10.8d0,10.8d0,10.8d0,10.8d0,10.8d0,10.8d0,16.99d0,
C          32      33      34      35      36      37      38    
     &   17.10d0,16.37d0,12.64d0,12.47d0,12.01d0,24.67d0,24.67d0,
C          39      40      41      42      43      44      45  
     &   24.67d0,24.67d0,24.67d0,24.67d0,24.67d0,24.67d0,24.67d0,
C          46      47      48      49      50      51      52
     &   24.67d0,24.67d0,24.67d0,37.32d0,38.71d0,38.44d0,31.74d0,
C          53      54      55        56         57
     &   31.50d0,29.99d0,315.275d0,226.994d0,176.252d0,
C          58      59       60       61        62       63       64
     &  140.68d0,140.68d0,140.68d0,140.68d0,140.68d0,140.68d0,140.68d0,
C          65      66       67       68        69       70       71
     &  140.68d0,140.68d0,140.68d0,140.68d0,140.68d0,140.68d0,140.68d0,
C          72
     &  105.112d0,
C          73      74       75    76      77      78     79
     &  81.24d0,81.24d0,81.24d0,81.24d0,81.24d0,81.24d0,81.24d0,
C          80      81       82       83       84       85       86
     &  57.364d0,57.254d0,63.162d0,63.540d0,55.283d0,57.171d0,56.64d0 /

      DO I = 1,N_ELEM_D2
        C600(I) = C6(I)
      ENDDO

      RETURN
      END


      SUBROUTINE SET_RCOV(N_ELEM,RCOV00)
      implicit none
      !Pass in as RCOV00 space for RCOV
      !Load from data statement 
      DOUBLE PRECISION RCOV(94)
      INTEGER N_ELEM,I
      DOUBLE PRECISION RCOV00(N_ELEM)
C
C     Covalent radii, original values Pyykko and Atsumi, Chem. Eur. J. 15, 2009, 188-197
C     Here scaled by k2 (=4/3) and converted to a_0 (factor = 0.52917726d0)
C     Elements H --> Pu indexed by Z.
C
      DATA RCOV/
C          1           2           3           4           5
     & 0.80628308, 1.15903197, 3.02356173, 2.36845659, 1.94011865,
C          6           7           8           9           10
     & 1.88972601, 1.78894056, 1.58736983, 1.61256616, 1.68815527,
C          11          12          13          14          15
     & 3.52748848, 3.14954334, 2.84718717, 2.62041997, 2.77159820,
C          16          17          18          19          20
     & 2.57002732, 2.49443835, 2.41884923, 4.43455700, 3.88023730,
C          21          22          23          24          25
     & 3.35111422, 3.07395437, 3.04875805, 2.77159820, 2.69600923,
C          26          27          28          29          30
     & 2.62041997, 2.51963467, 2.49443835, 2.54483100, 2.74640188,
C          31          32          33          34          35
     & 2.82199085, 2.74640188, 2.89757982, 2.77159820, 2.87238349,
C          36          37          38          39          40
     & 2.94797246, 4.76210950, 4.20778980, 3.70386304, 3.50229216,
C          41          42          43          44          45
     & 3.32591790, 3.12434702, 2.89757982, 2.84718717, 2.84718717,
C          46          47          48          49          50
     & 2.72120556, 2.89757982, 3.09915070, 3.22513231, 3.17473967,
C          51          52          53          54          55
     & 3.17473967, 3.09915070, 3.32591790, 3.30072128, 5.26603625,
C          56          57          58          59          60
     & 4.43455700, 4.08180818, 3.70386304, 3.98102289, 3.95582657,
C          61          62          63          64          65
     & 3.93062995, 3.90543362, 3.80464833, 3.82984466, 3.80464833,
C          66          67          68          69          70
     & 3.77945201, 3.75425569, 3.75425569, 3.72905937, 3.85504098,
C          71          72          73          74          75
     & 3.67866672, 3.45189952, 3.30072128, 3.09915070, 2.97316878,
C          76          77          78          79          80
     & 2.92277614, 2.79679452, 2.82199085, 2.84718717, 3.32591790,
C          81          82          83          84          85
     & 3.27552496, 3.27552496, 3.42670319, 3.30072128, 3.47709584,
C          86          87          88          89          90
     & 3.57788113, 5.06446567, 4.56053862, 4.20778980, 3.98102289,
C          91          92          93          94           
     & 3.82984466, 3.85504098, 3.88023730, 3.90543362 /

      DO I = 1, N_ELEM
         RCOV00(I)=RCOV(I)
      ENDDO

      RETURN
      END


      SUBROUTINE SET_R2R4(N_ELEM,R2R400)
      implicit none
      !Pass in as R2R400 space for R2R4
      !Load from data statement 
      DOUBLE PRECISION R2R4(94)
      INTEGER N_ELEM,I
      DOUBLE PRECISION R2R400(N_ELEM)
      DATA R2R4 /
C          1           2           3           4           5
     & 2.00734898,  1.56637132,  5.01986934,  3.85379032,  3.64446594,
C          6           7           8           9           10
     & 3.10492822,  2.71175247,  2.59361680,  2.38825250,  2.21522516,
C          11          12          13          14          15
     & 6.58585536,  5.46295967,  5.65216669,  4.88284902,  4.29727576,
C          16          17          18          19          20
     & 4.04108902,  3.72932356,  3.44677275,  7.97762753,  7.07623947,
C          21          22          23          24          25
     & 6.60844053,  6.28791364,  6.07728703,  5.54643096,  5.80491167,
C          26          27          28          29          30
     & 5.58415602,  5.41374528,  5.28497229,  5.22592821,  5.09817141,
C          31          32          33          34          35
     & 6.12149689,  5.54083734,  5.06696878,  4.87005108,  4.59089647,
C          36          37          38          39          40
     & 4.31176304,  9.55461698,  8.67396077,  7.97210197,  7.43439917,
C          41          42          43          44          45
     & 6.58711862,  6.19536215,  6.01517290,  5.81623410,  5.65710424,
C          46          47          48          49          50
     & 5.52640661,  5.44263305,  5.58285373,  7.02081898,  6.46815523,
C          51          52          53          54          55
     & 5.98089120,  5.81686657,  5.53321815,  5.25477007, 11.02204549,
C          56          57          58          59          60
     &10.15679528,  9.35167836,  9.06926079,  8.97241155,  8.90092807,
C          61          62          63          64          65
     & 8.85984840,  8.81736827,  8.79317710,  7.89969626,  8.80588454,
C          66          67          68          69          70
     & 8.42439218,  8.54289262,  8.47583370,  8.45090888,  8.47339339,
C          71          72          73          74          75
     & 7.83525634,  8.20702843,  7.70559063,  7.32755997,  7.03887381,
C          76          77          78          79          80
     & 6.68978720,  6.05450052,  5.88752022,  5.70661499,  5.78450695,
C          81          82          83          84          85
     & 7.79780729,  7.26443867,  6.78151984,  6.67883169,  6.39024318,
C          86          87          88          89          90
     & 6.09527958, 11.79156076, 11.10997644,  9.51377795,  8.67197068,
C          91          92          93          94           
     & 8.77140725,  8.65402716,  8.53923501,  8.85024712 /

      DO I = 1,N_ELEM
        R2R400(I)=R2R4(I)
      ENDDO

      RETURN
      END


      SUBROUTINE GET_COORD_CHAR(FCOORD,FCOORD2,CHARGE2)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
#include "pgroup.h"
#include "chrxyz.h"
#include "chrsgn.h"
#include "chrnos.h"

      DIMENSION FCOORD(NATOMS*3)
      DIMENSION FCOORD2(NATOMS,3)
      DIMENSION CHARGE2(NATOMS)

C   Unpack the coordinates and charges!
      IPCMN  = 0
      NCENTS = 0
C
      ICRX   = 1
      ICNT   = 1
      DO ICENT = 1, NUCIND
         MULCNT = ISTBNU(ICENT)
         IF (MULT(MULCNT) .EQ. 1) THEN
            IPCMN = IPCMN + 1
            NCENTS = NCENTS + 1
            FCOORD(ICRX  ) = CORD(1,ICENT)
            FCOORD(ICRX+1) = CORD(2,ICENT)
            FCOORD(ICRX+2) = CORD(3,ICENT)

            CHARGE2(ICNT) = CHARGE(ICENT)
            FCOORD2(ICNT,1) = CORD(1,ICENT)
            FCOORD2(ICNT,2) = CORD(2,ICENT)
            FCOORD2(ICNT,3) = CORD(3,ICENT)

            ICRX = ICRX + 3
            ICNT = ICNT + 1
         ELSE
            JATOM = 0
            DO ISYMOP = 0, MAXOPR
               IF (IAND(ISYMOP,MULCNT) .EQ. 0) THEN
                  IPCMN = IPCMN + 1
                  JATOM  = JATOM + 1
                  NCENTS = NCENTS + 1
                  FCOORD(ICRX  ) = PT(IAND(ISYMAX(1,1),ISYMOP))
     &                        *CORD(1,ICENT)
                  FCOORD(ICRX+1) = PT(IAND(ISYMAX(2,1),ISYMOP))
     &                        *CORD(2,ICENT)
                  FCOORD(ICRX+2) = PT(IAND(ISYMAX(3,1),ISYMOP))
     &                        *CORD(3,ICENT)
                  CHARGE2(ICNT) = CHARGE(ICENT)
                  FCOORD2(ICNT,1) = FCOORD(ICRX)
                  FCOORD2(ICNT,2) = FCOORD(ICRX+1)
                  FCOORD2(ICNT,3) = FCOORD(ICRX+2)

                  ICRX = ICRX + 3
                  ICNT = ICNT + 1
               END IF
            ENDDO ! ISYMOP  
         END IF ! MULTIPLICITY OF CENTRE
      ENDDO ! ICENT 

      IPRLOC=0
      IF (IPRLOC.GT.1) THEN
        CALL HEADER('Centre Info.',-1)
        WRITE(LUPRI,'(A10,1X,A10,3(1X,A20))')
     &             'Atom No.', 'Charge', 'X','Y','Z'
        DO I =1,NATOMS
          WRITE(LUPRI,'(I10,1X,F10.5,3(1X,F20.12))')
     &                I,CHARGE2(I),FCOORD2(I,1),
     &                FCOORD2(I,2),FCOORD2(I,3)
        ENDDO
      ENDIF

      RETURN
      END


      SUBROUTINE PACK_SYM_GRAD(TEMPG,NSIGN,NBASE,WORK,LWORK)
#include "implicit.h"
#include "priunit.h"
#include "maxaqn.h"
#include "mxcent.h"
#include "maxorb.h"
#include "cbinuc.h"
#include "nuclei.h"
#include "symmet.h"
#include "energy.h"
#include "dftcom.h"
#include "dftd.h"
#include "pgroup.h"
#include "chrxyz.h"
#include "chrsgn.h"
#include "chrnos.h"

      DIMENSION WORK(LWORK)
      DIMENSION TEMPG(NATOMS*3)
      INTEGER NSIGN(8),NBASE(8)

C     ********************************
C     ***** Symmetry Coordinates *****
C     ********************************
C
      ISTRT = 1
      IPREAD =1
      NCOOR  = 3*NUCDEP
      IF (LWORK .LT. 2*NCOOR*NCOOR) CALL STOPIT('PRICAR','TRACOR',
     &                                          2*NCOOR*NCOOR,LWORK)
      CALL TRACOR(WORK(ISTRT),WORK(ISTRT + NCOOR*NCOOR),1,NCOOR,IPREAD)
      IF (MAXREP .GT. 0) THEN
         CALL HEADER('Symmetry Coordinates',1)
         WRITE (LUPRI,'(A,8I5)')
     &            '  Number of coordinates in each symmetry: ',
     &            (NCRREP(I,1),I=0,MAXREP)
         DO ISYM = 0, MAXREP
         IF (NCRREP(ISYM,1) .GT. 0) THEN
            WRITE (LUPRI,'(/A,2X,A3,A2,I2,A1/)')
     &          '  Symmetry',REP(ISYM),' (',ISYM+1,')'
            DO IATOM = 1, NUCIND
               DO ICOOR = 1, 3
                  ISCOOR = IPTCNT(3*(IATOM - 1) + ICOOR,ISYM,1)
                  IF (ISCOOR .GT. 0) THEN
                     NB = 0
                     DO I = 1, NCOOR
                        IADR = ISTRT + (I - 1)*NCOOR + ISCOOR - 1
                        NSGN = NINT(WORK(IADR))
                        IF (NSGN .NE. 0) THEN
                           NB = NB + 1
                           NBASE(NB) = I
                           NSIGN(NB) = NSGN
                        END IF
                     ENDDO !NCOOR 
                     IF (NB .EQ. 1) THEN
                        WRITE (LUPRI,'(I5,3X,A,2X,A,I5)') ISCOOR,
     &                     NAMEX(3*IATOM)(1:4), CHRXYZ(-ICOOR), NBASE(1)
                           TEMPG(ISCOOR) = GRADFTD(NBASE(1))
                     ELSE
                        WRITE (LUPRI,'(I5,3X,A,2X,A,3X,A,I3,'
     &                   //CHRNOS(NB-1)//'(2X,A,I5),A,I1)')
     &                   ISCOOR,NAMEX(3*IATOM)(1:4),CHRXYZ(-ICOOR),' [',
     &                   NBASE(1),(CHRSGN(NSIGN(I)),NBASE(I),I=2,NB),
     &                   ' ]/',NB
                         ESUM =GRADFTD(NBASE(1))
                         DO I =2,NB
                            ESUM = ESUM + NSIGN(I)*GRADFTD(NBASE(I))
                         ENDDO
                         TEMPG(ISCOOR) =  ESUM
                     END IF
                  END IF
               ENDDO !ICOOR 
             ENDDO !IATOM 
         END IF
         ENDDO !MAXREP
      END IF

      IPRLOC=0
      IF (MAXREP.GT.0) THEN
        IF (IPRLOC.GT.1) THEN
        CALL HEADER('Sym. DFT-D Gradient Contribution',-1)
        IOFF = 0
        DO IATOM = 1, NUCDEP
            WRITE (LUPRI, '(1X,A6,F17.10,2F24.10)')
     &             NAMDEP(IATOM), (TEMPG(IOFF+J),J=1,3)
                   IOFF = IOFF + 3
        END DO
        ENDIF

        DO I = 1, NATOMS*3
          GRADFTD(I) = TEMPG(I)
        ENDDO
!Print out the gradient when symmetry is used by call to PRIGRD.
!Careful this work is not overwritten between the use above and here...
!Input: Gradient in Symm Coord. Output: Print Cart. Grad.
        IF (IPRLOC.GT.1) THEN
          CALL HEADER('Output from PRIGRD',-1)
          CALL PRIGRD(GRADFTD,WORK(ISTRT),WORK(ISTRT + NCOOR*NCOOR))
        ENDIF
      ENDIF

      RETURN
      END
